A review of recent progress in thermoelectric materials through computational methods

Reducing our overwhelming dependence on fossil fuels requires groundbreaking innovations in increasing our efficiency in energy consumption for current technologies and moving towards renewable energy sources. Thermoelectric materials can help in achieving both goals. Moreover, because of recent advances in high-performance computing, researchers more increasingly rely on computational methods in discovering new thermoelectric materials with economically feasible performance. In this article, significant thermoelectric materials discovered through these computational methods are systematically reviewed. Furthermore, the primary computational tools that aid the design of the next-generation thermoelectric materials are introduced and discussed. These techniques include various levels of density functional theory, electronic transport simulations, and phonon calculations.


Introduction
More than 70% of the global energy consumption is lost during conversion as waste heat [1], released to the environment at a wide range of temperatures (300-1000 K). This wasted energy originates from numerous sources such as mining and smelting plants, factories, electricity generators, air-conditioning and lighting appliances, pipelines, vehicles' engines, brakes and exhausts, and even the human body.
Harvesting this waste heat is of paramount importance for improving energy consumption efficiency. One promising solution is the thermoelectric (TE) effect [2], which governs the direct conversion of heat to electricity in a TE material. Additionally, TE materials can further reduce the CO 2 emissions, and enable extensive deployment of battery-and gridfree devices, wearable electronics, the Internet of Things (IoT), and much more. The performance of a thermoelectric material at a given temperature T is quantified by its figure of merit (ZT). ZT is defined as ZT = S 2 Coeff T∕ , where S Coeff is the Seebeck coefficient, σ is the electrical conductivity, and κ is the thermal conductivity. κ is made of two contributions: the electronic thermal conductivity (κ e ) and the lattice or phonon thermal conductivity (κ L ). Maximizing ZT in thermoelectric materials requires three criteria satisfied: first, a large S Coeff , to obtain the highest electric potential difference for a given temperature gradient; second, a low κ, to delay the thermal equilibration between the heat source and the environment; last, a high σ, to generate the highest current for the thermo-voltage.
To construct a TE harvesting unit, one must interface an n-type TE material to a p-type TE material, forming an n-p TE junction. Then, by putting one side of the TE junction in contact with a hot surface (waste heat source) and the other side with a cold surface (ambient), electricity is generated on the cold side. Several couples of the TE junctions are usually coupled together, both serially and in parallel, forming a TE module, shown in Fig. 1, to increase the amount of harvested energy. Such a configuration is commonly achieved by cutting n-and p-type bars from sintered pellets [3], or preparing n-and p-type thin films epitaxially deposited on an inflexible single-crystal substrate [4]. The primary appeal of TE generators is the fact that they are capable of converting waste heat directly into electrical energy without any moving parts. As a result, research thrived and plenty of TE materials (alloys, oxides, and nitrides) were discovered. After nearly a century of research and development, however, the practical use of TE devices is still limited to Fig. 1 Structural design of wearable thermoelectric harvesting device: a schematic view; b close-up view of one thermoelectric module. Adapted with permission from Ref. [12]. Copyright 2018 IEEE niche applications [5][6][7]. There are three significant obstacles against the widespread application of TE harvesters: (i) the low conversion efficiency and low chemical stability of TE materials; (ii) the toxicity and criticality of elements used in TE materials; and (iii) the mechanical rigidity of TE modules that impedes application in wearables.
In the search for high-performance thermoelectric materials, it has now become evident that the traditional materials' engineering approach based on judicious choice of materials and synthesis conditions via trial-and-error would only result in little improvement of the TE performance [8,9]. As a result, a theory-guided approach in the search for suitable high-performance thermoelectric materials with favorable electronic and thermal transport is becoming common [10][11][12]. In this review, we, therefore, systematically examine and present the recent developments in thermoelectric research propelled by theoretical and computational methods. The theoretical tools reviewed here can be divided into two broad categories: (a) electronic-structure calculations and that deal with S, σ, and κ e ; and (b) lattice dynamics or phonon calculations that determine κ L . Although these two types of calculations are interrelated, they, nonetheless, describe fundamentally different aspects of the materials' properties. As such, different computational tools have been developed for electronic and lattice dynamics calculations, of which both are covered here.

The need for computational insight
The second law of thermodynamics fundamentally dictates the upper limit efficiency of any thermodynamic engine, including the thermoelectric power generators. The efficiency (η) of a solid-state power generator, as a function of ZT, operating between hot and cold reservoirs with temperatures T H and T C is given by the following nonlinear expression [13]: The challenge in achieving high ZT ( S 2 Coeff T∕ ) materials lies in the interdependency of the κ, σ, and S Coeff in any given material system, as all are functions of the carrier concentration (n). For instance, in an isotropic n-type semiconductor, σ can be expressed as: in which e is the unit charge and μ e (T,n) is the electron mobility. μ e (T,n) is governed by a multitude of scattering mechanisms that often vary with T and carrier concentration, in a complex manner. n(T), in turn, is commonly controlled (2) (T, n) = e ⋅ n(T) ⋅ e (T, n), by extrinsic doping. S Coeff can be generally expressed by the formula [14]: in which k B is the Boltzmann constant, and N D is the concentration of n-type (donor) dopants. One should note that Eq. 3 only holds at relatively large n(T) as S Coeff diverges when n(T) approaches zero [15]. As can be easily deciphered, σ in thermoelectric semiconductors can be raised by carrier doping. As an undesired side effect, however, the additional carriers decrease S Coeff , resulting in a meager or no net benefit to the ZT.
Furthermore, the relation between κ and n is even more complicated. κ e is related to σ by the Wiedemann-Franz law [16]. As a result, any carrier doping unfavorably boosts the electronic share of the thermal conductivity, while κ L , which is the dominant thermal conduction mechanism over a wide range of carrier concentrations in semiconductors, depends on the density, heat capacity, and thermal diffusivity of which all are influenced by doping. For all these complications, κ is often phenomenologically modeled as a temperature-and concentration-dependent power law [17]: Here, κ 300 is the thermal conductivity at 300 K, while n ref , σ κ , and β κ are constants to be determined by fitting experimental measurements. In the intrinsic range, both electrons and holes contribute toward κ, constituting a bipolar contribution towards κ. The heavier an n-type semiconductor is doped, however, the more generous the κ e contribution is towards κ.
A reflection upon the complex interdependence of S Coeff , σ, and κ, one recognizes the challenges of designing highefficiency thermoelectric materials. The challenges arising from carrier doping are schematically presented in Fig. 2. In addition to the electronic requirements, one must also consider other constraints such as the non-toxicity and the abundance of the elements used (usually associated with the price and sustainability), mechanical robustness, moist resistance, flexibility (for wearables), and the low dopant solubility in some cases, to mention a few. As a result, the typical experimental approach in finding new thermoelectric materials inspired by search within a composition or structural space similar to known materials is very restrictive by nature. Consequently, by resorting to computationally guided materials search, design, and discovery, one may accelerate the pace of the progress in TE research beyond the limit imposed by the judicious and serendipitous experiments.

The electronic description
The efficiency of the thermoelectric energy conversion is directly related to the materials' non-dimensional figure of merit, which grows linearly with the electrical conductivity [21]. The essential quantities are then σ, together with the S Coeff , which are determined by the electronic structures of the materials and by carrier-scattering processes. Currently, there are two main approaches for calculating the electron transport properties and S Coeff : the semiclassical formulas based on Boltzmann transport theory; and the fully quantum-mechanical approach described in the Kubo-Greenwood theory [22][23][24]. An alternative, fully quantum-mechanical approach has also been developed by Landauer [25]. Here, the scattering processes are modeled as if they build charges and fields inside the materials, which represents an alternative view to that of Kubo's method in which currents are formed in response to a given external electric field. The Landauer methodology intimately resembles the experimental viewpoint, where one usually imposes an external current and measures the resulting potential drop due to the electron scattering, which is an important tool in guiding intuition when studying mesoscopic transport properties. The semiclassical Boltzmann transport equation, however, has gained appreciable popularity in recent years due to the ease of use and its integration with a diverse range of density functional theory (DFT) codes. The semiclassical Boltzmann transport equation, as implemented in the prominent BoltzTrap [24] and BoltzTrap2 [26] packages, relies only on the DFT calculated band and k-dependent quasiparticle energies (ε) as input. It then extrapolates the transport properties within the rigid band approximation [24]. Furthermore, this method assumes a T and ε independent relaxation time (τ) for electronic scatterings which results in a tractable and straightforward form of the equations for S Coeff , σ, and κ e . The constant relaxation-time approximation, despite its simplicity, predicts S Coeff values that match well with experiments and is widely adopted in the high-throughput theoretical search of novel TE materials. Table 1 summarizes some of the recent predictions using the BoltzTrap package.
In brief, the Boltzmann transport equation describes the general statistical behavior of a system out of thermodynamic equilibrium in terms of local equilibrium distributions. The Boltzmann transport equation is derived by assuming random motion of particles-in this case, electrons-within a medium, leading to the formulation of a current (j) in terms of conductivity tensors: in which σ ij and Ν ij are the elements of the electrical and thermoelectric conductivity tensors. From the band structure, within the rigid band approximation, the conductivity tensors along the reciprocal lattice direction α and β can be obtained as follows: where e is the electron charge, τ i,k is the relaxation time, v(i,k) is the group velocity, i is the band index, and k is the reciprocal (crystal momentum) vector. Here, v(i,k) can be calculated from the band structure by taking the derivative of the band energy function (f) over k. Here, the relaxation time (T i,k ) describes how quickly the system returns to thermodynamic equilibrium after scattering events. T i,k plays an important role in determining the transport properties such as electron mobility, electrical conductivity, thermal conductivity, and Seebeck coefficient. T i,k principally depends on both the band index (i) and the k vector direction. However, several studies of the direction dependence have shown that T i,k can be considered, with a good approximation, directionindependent in almost all known crystals [27]. A further approximation of T i,k to a constant value leads to an expression for the Seeback coefficient that is also independent of the relaxation time. However, an important restriction for the relaxation-time approximation to be valid is that the time the distribution function takes to relax to its equilibrium is independent of the distribution function itself and the externally

Fig. 2
The dependence of S Coeff , σ, and κ on carrier concentration (n) for polycrystalline Pb 1−x Sn x Te at 300 K. ZT reaches its maximum only for a very narrow window of n. A change in temperature shifts all these curves in a non-trivial manner, resulting in a maximum ZT at a different n. Furthermore, synthesis protocols may drastically change the presented values. S coeff and σ values were adapted from Refs. [18] and [19], while κ was adapted from Ref. [20] applied electric field. By integrating σ αβ over the energy (ε) space, the transport tensor can be obtained as: in which μ is the electronic chemical potential, and Ω is the unit cell volume. Finally, the Seebeck coefficient can be easily computed as follows: The physical limitation to the semiclassical Boltzmann theory is that it is only valid for the first-order terms, while it more or less fails at higher orders. This shortcoming is due to the assumption of chaotic molecular motion in Boltzmann's collision integral-an approximation that is only accurate up to first order [28]. However, the equivalence of the two most critical theoretical methods for the calculation of electronic transport within the linear regime has been explicitly demonstrated, and both of them allow the calculation of electric transport and S Coeff from first-principles (ab initio) calculations, namely DFT. The excellent agreement between the calculated results and the experimental data reported in the literature validates the validity of the transport theories when coupled with ab initio computational methods [29][30][31].

Wannier functions
The maximally localized Wannier function representation allows us to naturally introduce the ground-state electronic structure into the lattice Green's function approach that will be the basis for the evaluation of the Landauer quantum conductance [42,43]. This approach is computationally efficient, can be easily implemented as a postprocessing step in standard electronic-structure DFT calculations, and allows us to link the electronic transport properties of a device to the nature of the chemical bonds directly, providing an 300 S Coeff is maximum at T = 300 K, S Coeff is a convex function of T [41] insight into the mechanisms that govern electron flow at the nanoscale. Owing to the exponential localization of the Wannier functions obtained from planes waves, accurate integrations over the Brillouin zone can be achieved at a very moderate computational cost. Moreover, the analytical expression for the band derivatives in the Wannier basis resolves any issues that may occur when evaluating derivatives near band crossings. The accuracy of the results on a coarse k-point reciprocal-space grid using Wannier functions was confirmed by Pizzi et al., who performed calculations of the electrical conductivity on the binary and ternary skutterudites CoSb 3 and CoGe 3/2 S 3/2 [44]. Moreover, the semiclassical Boltzmann transport theory based on maximally localized Wannier functions with constant relaxationtime approximation was applied to examine the transport properties of SrTiO 3 [45]. The results for the transport properties show the prospect of thermoelectric presentation for SrTiO 3 at different temperatures from 300 to 800 K as a function of chemical potential, which plays an important role. Varying the temperature in terms of chemical potential shows a remarkable enhancement in the TE performance for the properties of thermoelectric materials. In the considered cubic phase, the chemical potential of SrTiO 3 is equal to the sum of all atomic component chemical potentials in the crystal [45].

In search of flat bands
The band dispersion, in other words, the slope of the bands near band edges, is of great importance for predicting materials' TE performance. Steeper slopes indicate stronger orbital interactions and faster carrier mobility. On one hand, for stronger orbital interactions, the band becomes more disperse in energy. On the other hand, weaker orbital interactions correspond to flatter bands. This simple model elucidates the fact that dispersed bands lead to higher carrier mobility, while flat bands indicate more localized carriers. However, there is a class of fermion systems with a diverging density of states (DOS). These type of materials, such as CuAlO 2 [46], exhibit a dispersionless spectrum that have precisely zero dispersion, the so-called flat band [47]. Flat bands can emerge in strongly interacting condensed matter systems and layered compounds with integer-valued pseudospin [48][49][50][51][52]. Ishii et al. found that the electric current flowing along the chain suddenly increases the magnitude when a minimal electric field is applied to the system. This large current originates from the symmetry breaking of the large energy degeneracy of flat-band states in the Kagome lattice, which is nearly independent of the size of the Kagome-lattice chain itself [53]. Under flat-band conditions, Ishii et al. showed that carriers in high-lying minibands could propagate for several superlattice periods without being scattered. These carriers follow the miniband description of transport. In a uniform electric field, the minibands break down into localized states immediately. These states participate in conduction. One strategy to achieve large S Coeff and large electron transport is via low-energy electronic bands containing both, flat and dispersive parts, in different regions of crystal momentum space. This approach is known as the pudding-mold band structure, of which examples are shown in FIG. 3 [54,55].
When we consider the problem of electron transport from a more general viewpoint, the literature suggests that most layered compounds are particularly interesting because of their unusual electronic and physical properties, in addition to their desirable thermoelectric performance [56]. For example, through electronic-structure calculations using DFT, a large and anomalous anisotropic thermoelectric transport properties were predicted in d 0 -electron layered complex nitrides AMN 2 (where A and M are metal cations). Subsequently, it was proposed that these compounds may constitute a novel and promising class of TE materials [57][58][59][60]. On one hand, because of the dimensionality of the electronic structures, anisotropic electronic transport properties have also been observed in layered oxides La x Sr 1−x CuO 4 and Sr 2 RuO 4 [61][62][63]. On the other hand, complex layered nitrides in the form AMN 2 (A = Sr or Na; M = Zr, Hf, Nb, Ta), were analyzed. The Sr containing compounds SrZrN 2 and SrHfN 2 were found to have three-dimensional electronic structures and isotropic thermoelectric transport properties, despite their α-NaFeO 2 -type layered crystal structure [64]. This apparent contradiction can be explained by the fact that the bottom of the conduction bands in these compounds is composed of Zr 4d z 2 + 4d xz + 4d yz and Hf 5d z 2 + 5d xz + 5d yz , along with Sr 4d x 2 −y 2 + 4d xy orbitals. Here, the A and M ions jointly provide a full set of the d orbitals. However, although the electronic structures of NaNbN 2 and NaTaN 2 are also three-dimensional, NaNbN 2 and NaTaN 2 were found to have weak anisotropic thermoelectric transport properties. This weak anisotropy arises from the NbN 2 and TaN 2 layers whose conduction band minimums are composed of Nb 4d z 2 + 4d xz + 4d yz and Ta 5d z 2 + 5d xz + 5d yz along with N 2p z orbitals with no contribution from the Na orbitals. Understanding and elucidating the origin of anisotropic transport properties in layered compounds are of great significance for cultivating this new class of TE compounds, not only for thermoelectric materials but also for realizing novel physical properties [65].

Spin-driven Seebeck effect
Most conventional thermoelectric materials, of which we reviewed some in the previous section, are semiconductors, rather than metals, and, consequently, have low intrinsic electrical conductivity. That is because although metals are the ideal conductors, they, nonetheless, have a diminishingly small S Coeff . The carrier symmetry at the Fermi level (E F ) can explain metals' small Seebeck coefficient, as it implies the number of hot electrons above E F and cold empty states below E F are roughly equal. When the number of diffused holes and electrons under a temperature difference is roughly the same, the net Seebeck effect is minimal due to cancelation. Semiconductors circumvent this problem by having only one type of majority carriers. There is, however, one class of thermoelectric materials that noticeably deviates from the conventional paradigm, i.e., the strongly correlated thermoelectrics. Thermoelectric materials such as Ca 3 Co 4 O 9+δ [66], Fe 1+y Te 0.6 Se 0.4 [67], and Na x CoO 2 [68] are excellent electrical conductors-Na x CoO 2 has a high carrier concentrations (n) which is in the order of ~ 10 21 to 10 22 cm −3 [69]. In these materials, spin entropy flow drives the Seebeck effect. For instance, in Na x CoO 2 , the susceptible spin population is approximately equal to the hole population and remains constant with decreasing T. Na x CoO 2 , therefore, has a conductivity similar to a metal and a magnetic behavior similar to a frustrated insulator. Such compounds are known as Curie-Weiss metals [70,71]. Wang et al. [72] demonstrated that at T = 2 K, a magnetic field completely suppresses S Coeff to nil. The disappearance of S Coeff in Na x CoO 2 , when an external magnetic field aligns all the spins, implies that the Seebeck effect originates from the spin entropy flow carried by the holes in the Curie-Weiss phase. In this section, therefore, we review the structural and electronic characteristics of Na x CoO 2 -the most Fig. 3 The red curves in a, b show pudding-mold band structure in CuAlO 2 [46] and BaPdS 2 [54], respectively. Schematic figure of the electronic structure within the pudding-mold structure shown in (c) and compare to the one of a typical metal (d). Moreover, e, f represent the band dispersion [55]. Here, t 1 is the hopping integral between neighboring sites, and ε is the quasiparticle energy. The dotted gray lines mark the position of the chemical potential for a given temperature T. Adapted with permission from Refs. [46,54,55]. Copyrights 2013, and 2019 American Physical Society and 2007 Physical Society of Japan investigated such compound-as a representative of the strongly correlated thermoelectric material.

Na x CoO 2
Sodium cobaltate (Na x CoO 2 ) is a strongly correlated oxide [73,74] with a fairly high thermoelectric figure of merit (ZT) at temperatures around 800 K for x > 0.5 [75]. As shown in Fig. 4, the Na x CoO 2 lattice is made of alternating Na and edge-sharing CoO 2 octahedral layers stacked along the c direction constituting a quasi-two-dimensional system. In Na deficient (x < 1) Na x CoO 2 , each Na vacancy creates a hole converting a Co 3+ ion into Co 4+ . As a result, Na x CoO 2 exhibits p-type conductivity with n inversely proportional to x. Moreover, the coexistence of Co 3+ and Co 4+ ions, or charge disproportionation, in the CoO 2 layers generates a competitive Seebeck potential through spin entropy flow [68,72]. Furthermore, the mobile Na + ions diffusively scatter the phonons [70,76], decreasing κ L to ~ 0.01 W cm −1 K −1 at ~ 1000 K [77]-this κ L value is ~ 100 times smaller than that of ZnO at the same temperature [78].
The electrical transport properties of Na x CoO 2 have also been reported for different concentrations of Na in several experimental studies [70,79,80]; however, the results are inconsistent at times and even contradictory. For example, Kawata et al. [79] found that for x > 0.55, the electrical resistivity (ρ) of Na x CoO 2 increases along with increasing x. Motohashi et al. [80], however, found the opposite trend for ρ as a function of x. More recently, Foo et al. used the chemical de-intercalation from Na 0.77 CoO 2 to change the Na concentration and found a metal-insulator transition of the charge ordering in Na 0.5 CoO 2 single crystal. A similar effect, however, could not be detected in the polycrystalline Na x CoO 2 prepared by solid-state reaction. In another work, Cui et al. [81] reported a comprehensive study of the structure and electrical transport properties of Na x CoO 2 (x = 0.52, 0.56, 0.60, 0.68, 0.72, 0.76, 0.8). Here, all the samples had metallic conduction, with ρ reaching a minimum at around x = 0.64. Moreover, a jump in ∂ρ/∂T was observed at the temperature corresponding to the dimensional crossover of the electronic transport properties. Some of the discrepancies can be explained by the structural differences between Na x CoO 2 prepared by solid-state reaction and by chemical de-intercalation from Na 0.77 CoO 2 . For instance, it was shown that the distance between Co in the CoO 2 layers, which may depend on the synthesis protocol, is critical in determining the electrical transport property of Na x CoO 2 -and more generally, of all the Na x TMO 2 (TM = Cr, Mn, Fe, Co, Ni)where this critical distance R c that determines whether the compound is insulating or metallic [82]. As seen here, the origin of this difference in observations is still an open question [83,84] and sets a challenge for theoretical analysis.

Spin entropy flow
The high-temperature S Coeff in Na x CoO 2 can be expressed by Koshibae's equation [85,86] for strongly correlated materials, which itself is a modified form of Heikes formula [87,88]: Here, g refers to the degeneracy of the electronic configuration in ions A and B (here A = Co 3+ and B = Co 4+ ), and it is equal to all the distinct possible ways which the electrons can be arranged in the orbitals of the ions they occupy. x is the ratio of Co 4+ over all Co ions, which is proportional to the carrier (hole) concentration, and finally, k B and e are the Boltzmann's constant and the electron's charge, respectively. This concept is schematically demonstrated in Fig. 5. In Na x CoO 2 , the six electrons of the low-spin Co 3+ have a t 6 2g e 0 g configuration for which there is no other degenerate state that is entropy S = ln[g(Co 3+ )] = ln1. In contrast, for the lowspin Co 4+ ion, since one electron is removed, there are six d e g e n e r a t e c o n f i g u r a t i o n s . C o n s e q u e n t l y, S = ln[g(Co 4+ )] = ln6. As an electron hops from a Co 3+ ion to a Co 4+ ion, an entropy flux moves opposite to the electric current flow (hole); this phenomenon is referred to as spin entropy flow [89]. According to Eq. 11, the S Coeff for Na 0.75 CoO 2 (a common representative of the Na x CoO 2 family with metallic conduction) is 249 μV K −1 , which is quite comparable with the experimental measurements of ~ 200 μV K −1 at 800 K [90]. Xiang et al. attempted to explain Na x CoO 2 's high S Coeff value using standard Boltzmann transport theory combined with spin-polarized DFT [92]. They found that S Coeff was smaller when the system was set to be magnetically polarized, and interpreted the results as an alternative explanation for the suppression of S Coeff in a magnetic field. According to their interpretation, the high S Coeff and its suppression in high magnetic fields could be understood using the conventional transport theory applicable to other thermoelectric materials and that the unique thermoelectric properties of Na x CoO 2 are due to its unusual band structure as manifested in the narrow manifold of t 2g bands in this system. However, the dependence of S Coeff on Co's spin state (high spin vs. low spin) [93] still points to the spin entropy flow as the primary cause for high S Coeff values.

What level of DFT?
Density functional theory (DFT) is one of the most common theoretical tools for probing strongly correlated TE materials such as Na x CoO 2 . Within the DFT framework, the most basic approximation for the exchange-correlation functional, i.e., the local density approximation (LDA), in which the exchange and correlation energy terms are fitted to those of a homogenous electron gas [94], fails to describe the electronic structures in highly correlated systems realistically [95]. The next most popular functional based on the general gradient approximation (GGA) [96], often, offers only a minor improvement [97]. The shortcoming stems from the fact the both LDA and GGA functionals tend to delocalize electrons over the lattice space. As such, each electron feels an average of the Coulombic potential. Consequently, for highly correlated materials, such as Na 0.75 CoO 2 , the significant Coulombic repulsion between localized electrons that results in charge disproportionation, are not correctly represented by LDA and GGA [83]. The thermoelectric effect, nonetheless, stems from charge disproportionation and needs an accurate theoretical description.
A standard solution for avoiding the delocalization problem of LDA and GGA functionals is to add an ad hoc Hubbard-like localized term to the LDA and GGA functionals. This approach is known as DFT + U (commonly referred to as LDA + U and GGA + U). Practically, in the DFT + U method [98,99], one applies a Hubbard (U) and an exchange (J) terms to the localized d (or f) electrons. The delocalized electrons (s and p) are left unaltered as LDA and GGA correctly describe these delocalized electrons. To avoid the double-counting of the correlation part for localized electrons, another term-called the double-counting correction-is subtracted from the Hamiltonian. U − J, known as U eff , determines the level of localization within the DFT + U implementation.
The U eff value is not usually known a priori. Therefore, it is often chosen to replicate some experimental measurements [100]. For this reason, multiple values are usually reported in the literature for a given compound depending on the experimental values used as a benchmark, such as the band gap, formation enthalpy, or magnetic phase. In the case of Na x CoO 2 , among the reported U and J values in the literature [101][102][103], all reproduce the charge disproportionation [104,105] of the Co ions more or less accurately. That is to clearly distinguish Co 2+ , Co 3+ , and Co 4+ from one another in terms of magnetization and partial density of states (PDOS). For instance, the more common value of U eff = 4 eV is at a midpoint between the lower values of ~ 3 eV implemented Co 3+ [106] and the higher values of ~ 5 eV implemented Co 4+ [107][108][109].

The criticality of U eff
As presented in the previous section, the exact choice of U eff for Na x CoO 2 remains somehow ambiguous. Here, we examine how the results of the DFT calculations with various U eff value for Na 0.75 CoO 2 (as a representative for the Na x CoO 2 compound family) differ from one another. The problem to address is the arrangement of Reproduced with permission from Ref. [91]. Copyright 2020 Institute of Physics Na ions (or equally Na vacancies) in the unit cell. Figure 6a shows the Na 0.75 CoO 2 obtained with the LDA level of the theory [110]. In this structure, the unit cell contains a √ 3a × 4a × 1c supercell in which the Na ions were arranged in rows of 3 × Na2, progressively followed by a Na vacancy, a 3 × Na1, and another Na vacancy. This structure has equal numbers of Na2 and Na1 ions (Na1/ Na2 = 1). Figure 6(b) shows the Na arrangement obtained with U eff = 4 eV (GGA + U formalism) [111]. Here, the irreducible unit cell has dimensions of 2 √ 3a × 2a × 1c in which one-third of the Na ions occupied Na1 sites while two-third occupied Na2 sites (Na1/Na2 = 0.5) as the Na vacancies create a multi-vacancy cluster as marked by green enclosures. Moreover, Na ions at the Na2 sites were ever slightly off the oxygen axes conforming to the H1 structure [112]. Figure 6c shows the Na arrangement in Na 0.75 CoO 2 calculated with U eff = 5 eV (GGA + U formalism) [113]. Here, the ratio of Na1/Na2 is ~ 0.44, where the Na ions in Na1 sites form ordered tri-Na1 droplets. The margin of the stability of the structure at Fig. 6c was, however, only ~ 2 meV lower than the competing structures which the authors had considered, implying that phase separation might occur at room temperature.
We can infer that the choice of U eff is demonstrably critical in describing the arrangement of Na ions in Na 0.75 CoO 2 . Consequently, to validate the theoretical results, an accurate comparison with the experiment is required. This task is, however, easier said than done. Na is volatile, and its concentration in the sample may be smaller than the nominal value [114]. Furthermore, deciphering the real atomic environment at a given lattice point in Na x CoO 2 experimentally, which is strongly dominated by the coordination environment, can be determined using the analysis of X-ray absorption near-edge spectra. However, such a measurement often is limited by the energy range of the instruments, rendering the interpretation of the results less than deterministic [115]. Obtaining a well-calibrated U eff value for Na x CoO 2 , therefore, remains a difficult task. Utilizing higher level theory such as Hartree-Fock/DFT hybrid functionals [116] or GW approximation [117] may offer more robust theoretical tools for probing Na x CoO 2 and similar materials. Unfortunately, these methods are more computationally demanding than the standard or + U corrected DFT calculations.

Phonon calculations
As mentioned above, the efficiency of a TE generator depends on the ZT of the TE material, where κ can be further separated into lattice (κ L ) and electronic (κ e ) thermal conductivities. The latter is correlated with σ as shown by the Wiedemann-Franz law: where L is the Lorenz number. Since κ L is the only independent variable in ZT, lowering κ L has been the most effective way to increase ZT in the past decades. An ideal TE material has the so-called "phonon-glass electroncrystal" (PGEC) property [118], with a high power factor (12) e ∕ = LT, Fig. 6 The schematic presentation of Na ion arrangements in Na 0.75 CoO 2 calculated with different U eff values. a LDA calculations based on Zhang et al. [110]. The open circles denote those Na ions at Z = 0, while solid circles denote those Na ions at Z = 0.5 within the unit cell. b GGA + U (U eff = 4 eV) calculations based on Assadi et al. [111]. Here, the color scheme is similar to that in Fig. 4. c GGA + U calculation based on Hinuma et al. [113] with U eff = 5 eV. Here, the Na ions at Z = 0 are at mirror locations of those at Z = 0.5. The blue and gold spheres rather denote Na1 and Na2. Adapted after Refs. [110,111,113]. Copyrights 2005 American Physical Society, 2015 World Scientific, and 2008 American Institute of Physics ( PF = S 2 Coeff ) and a low κ L . In this section, we will briefly summarize recent advances and challenges related to the ab initio description of κ L .

Phonons and thermal conductivity
Debye first used quantized normal mode of the collective atomic vibrations to describe the specific heat of crystalline solids, giving birth to the concept of the phonon [119]. The phonon description of lattice vibrations has been spectacularly successful and has become a pillar of solid-state physics. In a harmonic oscillating lattice, phonons never interact, there is no lattice expansion, and κ L approaches infinity. In a real lattice, the interatomic potential has anharmonic terms, which leads to interactions between phonon quasiparticles, thus resulting in a finite κ L and lattice thermal expansion (some intriguing materials exhibit negative thermal expansion within a temperature range [120][121][122][123]). Peierls adapted the Boltzmann theory for the transport of gases to phonon "gas" to describe their population distributions driven out of equilibrium by a temperature gradient, formulating the Peierls-Boltzmann transport (PBTE) equation for determining phonon lifetimes and κ L [124,125]. The solution of the PBTE, a complex set of coupled integrodifferential equations, is a complicated task by itself. A variety of approximations have thus been introduced to develop an understanding of phonon driven thermal properties. For example, one can mention the Callaway and Klemens models of intrinsic phonon scattering and conductivity [126,127]. Furthermore, Slack et al. developed a simple model for κ L and a simple set of rules for its interpretation based on crystal complexity, average atomic mass, bonding strength, and anharmonicity [128,129]. These are quantified in terms of a material's Debye temperature and average Grüneisen parameter, which gives a measure of the speed of sound of the heat-carrying phonons and anharmonicity, respectively.
With the recent advances in theoretical and computational techniques and the exponentially increasing computational power, phonon band structure, group velocities, lifetimes, and thermal conductivity can be calculated from ab initio without ad hoc adjustable parameters. The solution of the PBTE has been achieved so substantially beyond the typical relaxation-time approximation, from variational methods [130] to self-consistent iterative procedures [131] and compressive sensing [132]. Moreover, realistic interatomic potentials for a variety of materials can be obtained from DFT [133]. The harmonic and anharmonic interatomic force constants (IFCs), as in Fig. 7 and Eq. 13, can be determined numerically either from DFT electronic-structure calculations based on supercell perturbations method [134] or from the linear response methods within density functional perturbation theory (DFPT) [135]. In the anharmonic case, as shown in Fig. 7, the interatomic force is not proportional to the displacement, due to the third, fourth, and even higher order terms in the Taylor expansion of the potential energy V: Nowadays, fully ab initio calculations of lattice thermal conductivity have become standard practice, owing to the availability of computational power and open-source code packages including but not limited to ShengBTE [136], Pho-no3py [137], PhonTS [131], and almaBTE [138]. However, these standard computations represent perfect and infinite single crystals. Real materials are finite and with imperfections, point defects like vacancies and interstitials, dislocations, and grain boundaries. In TE materials, defects are typically desired to tune the electronic band structure via doping and provide important scattering centers for phonons to lower κ L . Defects add considerable complexity to the phonon calculations [139] and the prediction of κ L in realistic materials. Ab initio calculations of phonon-phonon and phonon-defect interactions on the same footing are just becoming possible recently.

Realistic phonon calculations
A single site point defect changes the local structure and chemical bonding around it, providing a scattering center for phonons decreasing the thermal conductivity. This Fig. 7 Comparison of a harmonic and an anharmonic oscillator thermal resistance mechanism can be accounted, to some extent, within the above-described software packages for calculating κ L , often demonstrating excellent agreement with the measured data for high-quality crystals. However, the theory assumes that the point defects are randomly distributed in the lattice within the diluted limit, meaning that defects can scatter phonons without any correlation between them. As seen in Fig. 8, in the engineered materials with high point-defect concentration, this independent scattering model is questionable [140][141][142]. The inadequacy of this first-order perturbative description of point-defect scattering was shown in the recent work when comparing results from the first-order approximation with those of fully ab initio Green's function calculations of phonon scattering by vacancies and nitrogen substitution in diamond [141]. This work showed that including the variation of IFC induced by point defects and going beyond the first-order perturbations give much more reliable phonon-defect scattering prediction than the simple mass perturbation model.
Further works have combined first-principles Green's function and PBTE methods to describe and predict κ L in a variety of materials with different defect types (intrinsic vacancies and antisite, and substitutional dopants) [143,144]. Among the intrinsic defects, vacancies are generally found to have the most potent effects on phonons scattering, due to the magnitude of the missing mass from vacancy centers and structural distortions. However, the "dilute limit" approximation is still assumed in most of the typical stateof-the-art phonon-defect calculations. Various methods have been employed for high concentrations of defects and coherent scattering mechanisms, such as the virtual crystal approximation (VCA) [145], supercell unfolding [142], and ab initio molecular dynamics (AIMD) simulations [146,147].
One avenue in the search for low κ L and high ZT is to examine very complex materials [148]; the ones that possess a variety of bonding patterns, rattling atoms in covalent cage structures, and many low-frequency optic phonon branches for the scattering of heat-carrying acoustic phonons. Computationally, such crystal complexity brings about significant challenges. In fact, ab initio PBTE methods have been recently used to explore complex unit cells of, for example, PbRb 4 Br 6 or Cu 12 Sb 4 S 13 [132,149]. The numerous vibrational degrees of freedom in the complex unit cells make the DFT calculations of the IFCs computationally tricky and increase the phase space for phonon-phonon interactions. These numerical bottlenecks would become less problematic as the computational power increases, offering a better understanding of the relevant physical mechanisms governing κ L , and allowing better approximations to be made without much loss of precision in the near future.
A part of the state-of-the-art ab initio thermal transport methodologies is based on the phonon quasiparticle scattering picture within the quantum perturbation theory. Intrinsic phonon-phonon interactions are derived from the perturbations coming from higher order anharmonicity of the interatomic potential (derivatives of the interatomic potential with respect to each atomic displacement, see Eq. 13). For materials with extreme behaviors (such as strongly anharmonic) or in extreme conditions [like high T or pressure (P)], the ab initio methods derived from lowest order perturbation theory could show discrepancies between measured and calculated phonon transport properties. Materials that are about to or have undergone phase transition around or below the temperature of interest pose a crucial challenge as their lattice is highly anharmonic. Explicitly incorporating the thermal effects onto the anharmonic IFCs, beyond the quasiharmonic approximation, therefore, becomes necessary. Furthermore, including the fourth-order anharmonicity has been shown to be essential for κ L in Si [150] at high temperature, and in PbTe, even at room temperature [151]. Molecular dynamics (MD) simulations of phonon transport include all orders of the anharmonicity. MD, however, treats the atomic interactions classically; and AIMD remains prohibitively costly to tackle realistic system sizes.
Describing the lattice vibration and thermal transport at extreme T and P and near phase transitions driven by T or P will be the challenge for future ab initio phonon transport calculation and even to DFT in general. Phonons in materials near phase transitions can have localization as the mean-free path becomes comparable to their wavelength, constituting the Ioffe-Regel regime [152]. In this circumstance, a wave can no longer be readily defined. A localized vibrational state can become hybridized with propagating phonon states, resulting in resonant-like extended states with diffusive motion, not propagating, in nature. Perhaps, new quasiparticle picture like diffusons/locons [153,154] in disordered materials will be needed to describe the lattice vibration phenomena in these extreme cases.
Recently, Hellman et al. have developed an accurate method, called the temperature-dependent effective potential (TDEP) technique, to determine the temperature-dependent anharmonic free energy, based on ab initio molecular dynamics followed by a mapping onto an effective model Hamiltonian describing the lattice dynamics. This method has been applied to various materials with dynamic instability, such as bcc Zr [155], and some strongly anharmonic materials, such as PbTe [156] and SnSe [157], which happen to be good thermoelectrics due to the very low thermal conductivity. Moreover, Zhou et al. calculated the electron-phonon interaction in the presence of soft-modes in SrTiO 3 perovskite from ab initio using (anharmonic) TDEP phonons [158]. The stochastic self-consistent harmonic approximation (SSCHA) is another novel development in the theoretical framework to study from ab initio the anharmonic properties of solids. SSCHA has been successfully applied to study the second-order structural phase transition and temperature-dependent anharmonic phonons primarily in various materials, such as the ferroelectric SnTe and GeTe [159], and the thermoelectric SnSe [159], as shown in Fig. 9.
For advanced TE materials, high-performance improvement needs to combine optimization of both the electronic and vibrational properties of the materials as these are tightly coupled. Thus, state-of-the-art numerical simulations of the TE materials need to incorporate electron-phonon interactions to improve our understanding of κ L S Coeff , and σ as coupled properties. Numerical techniques have been built and developed recently for efficient calculation of electron-phonon interaction matrix elements from DFPT and quantum perturbation theory [160,161]. These ab initio electron-phonon coupling has been examined and used to study transport phenomena in doped TE semiconductors [158,[162][163][164][165] and metals [166][167][168]. Rather than solving the transport equation for electron and the phonon separately, the solution of their coupled distribution equations would be the focus of future calculations. These calculations shall be critically important in understanding the phonon-drag effects in TE materials, and the highly out-of-equilibrium electron and phonon dynamics in the ultrafast laser-based thermal spectroscopies.

Interfacial systems
Effective control of interfaces in thermoelectric materials is critical for improving their efficiency. In general, interfaces in thermoelectric devices are divided into two main types: (a) grain boundaries within single-phase materials and solid solutions; and (b) interfacial systems formed between different materials, e.g., superlattices or epitaxial layers. Figure 10 shows a schematic representation of grain boundaries in polycrystalline materials and multilayered heterostructures. Grain boundaries constitute the interfaces between grains in polycrystalline bulk materials. These grains may exhibit different atomic structures, crystal orientations, or different compositions. The role of grain boundaries is especially crucial in nanostructured materials, which are composed of small crystal grain sizes and exhibit a high density of grain boundaries. However, despite the numerous experimental observations of grain boundaries, their implications on transport properties are not entirely understood. Theoretical simulations are, therefore, a powerful tool to study the role of grain boundaries and interfacial systems in thermoelectric materials at the nano-and sub-nano-scale. Computer simulations can be used to discuss the stability of specific interfacial systems, in which vacancies or defects are commonly present (and serve as scattering centers, as explained in the previous section), the materials' behavior under mechanical The calculations were performed with LDA on the experimental structure. Adapted with permission from Ref. [159]. Copyright 2019 American Society of Physics strain or shear deformations, and their corresponding TE performance. Such theoretical investigations are, however, computationally expensive and have just recently become popular within the thermoelectric community. In this section, we will review some of the landmark studies in the field, focusing on the computational techniques employed.

Nanostructuring and grain boundaries
From the theoretical perspective, electron and phonon transport in polycrystalline one-dimensional graphene nanoribbons and nanotubes was investigated by density functional tight-binding (DFTB) theory combined with non-equilibrium Green functions by Lehmann et al. [169]. The introduction of lattice imperfections such as dislocations was modeled by pentagons, heptagons, and octagons embedded within the hexagonal carbon nanostructures (Fig. 11a) which is suggested as an efficient method for enhancing the TE effect. The computed S Coeff exhibit values ~ − 0.4 mV/K and ~ − 0.6 mV/K in polycrystalline nanoribbons and nanotubes, respectively, in contrast to S Coeff > 0 in the ideal graphene structures. Variations of the electronic and thermal properties in polycrystalline carbon-based systems give rise to an improved ZT by several orders of magnitude at room temperature.
Phosphorene (composed of a few-layer-thick black phosphorous) is another bidimensional material that exhibits highly anisotropic electrical and thermal conductances in single-crystal form. A combined DFT and self-energy correction (GW) study [170] demonstrated that the preferred directions for κ and σ in a phosphorene monolayer are orthogonal to one another. This results in an anisotropic ZT  Structures of a polycrystalline graphene nanoribbon with pentagon (blue shaded area), heptagon (orange), and octagon (yellow) defects obtained by joining two layers with an inclination angle of 24° [169], b Bi 2 Te 3 one quintuple layer-thick slab (Te and Bi atoms are presented with orange and purple spheres, respectively) [173], and c p-type LaNiO 3 /SrTiO 3 (001) superlattice with two p-type interfaces (La, Sr, and O atoms are purple, green, and red, respectively. Light and dark blue octahedra are centered around Ti and Ni ions, respectively) [180]. Adapted with permission from Refs. [169,173,180]. Copyrights 2015 and 2017 American Physical Society, and 2010 American Institute of Physics that is large along the armchair direction, with ZT above 1 at 300 K and above 2 at 500 K, both for 60 ps relaxation time. However, in the experimentally processed material, the presence of defects may reduce the conductance and diminish the anisotropies.
In Si nanowires, simulations carried out by the ensemble Monte Carlo technique (EMC) [171] showed that the electrical conductivity decreases along with the cross section of the wire, while, on the other hand, S Coeff increases as the wire cross section decreases. The calculated ZT in nanowires is 20-40 times more substantial compared to the bulk Si. The decrease in the κ L due to strong phonon-boundary scattering can explain the remarkable enhancement in ZT. The abovedescribed studies show that the formation of grain boundaries can be used to improve the thermoelectric efficiency by maximizing S Coeff while decreasing the resistivity (ρ) and κ. On the other hand, grain boundaries can also affect the conductance negatively, reducing the anisotropy observed in single-crystalline counterparts.

Interfaces and heterostructures
Besides the recent progress in improving the performance of thermoelectric materials by modifying grain boundaries in nanostructured compounds, superlattices (formed by alternate stacking of two materials) have also shown potential to enhance the efficiency of TE energy conversion devices [172]. Ab initio DFT calculations showed that the slab formed by one atomic quintuple layer of Bi 2 Te 3 (Fig. 11b) has the capacity to increase ZT to 7.15, which is ten times larger than for the bulk structure at room temperature [173]. The significant enhancement in ZT is caused by the changes in the valence band density from quantum confinement within the thin film.
Saha et al. studied HfN/ScN [174] and ZrN/ScN [175] metal/semiconductor superlattices constructed by stacking layers of both materials alternatively. The structure of these model interfaces and the computed electrostatic potentials as a function of the distance are depicted in Fig. 12. The p-type Schottky barrier height for the HfN/ScN of 0.77 eV is larger than the 0.34 eV in ZrN/ScN superlattices. The n-type barrier of the HfN/ScN superlattices is estimated to be 0.13 eV (versus 0.56 eV for ZrN/ScN), which suggests that these superlattice models are degenerate for electron transport along the cross plane. The lower barrier in HfN/ ScN indicates the metallic conduction and lower S Coeff compared to the ZrN counterpart. On the other hand, HfN/ScN superlattices present better κ L than the ZrN-based models. These results point out the trade-off between S Coeff and σ, suggesting that simultaneously increasing both parameters is still challenging in two-dimensional geometries.
Ab initio calculations have also been used to discuss the stability of solid-state TE with high cooling power density, such as epitaxial Ni/Bi 2 Te 3 , NiTe/Bi 2 Te 3 , Co/Bi 2 Te 3, and CoTe 2 /Bi 2 Te 3 interfaces [176]. In these models, the metallic part is compressed/expanded by up to ± 10% from its equilibrium structure to match the Bi 2 Te 3 lattice parameters. All these interfaces were found to form stable Ohmic contacts. Other ab initio studies of epitaxial CoSb 3 /TiCoSb [177] and CoSb 3 /Ti [178] interfaces discussed the stability, electronic structure, and mechanical properties of these model structures. These publications [177,178] examined the atomiclevel mechanical strain failure mechanisms in CoSb 3 /Tibased interfacial systems and found that the fracture takes  [175] or c ScN and HfN [174]. b and c represent the planar average electrostatic potential (green dashed line) as a function of perpendicular distance from the ScN/ZrN and ScN/HfN interfaces, respectively. The red dashed line indicates the lattice-plane oscillations, which are filtered with the macroscopic averaging technique [182]. The vertical dashed black line represents the lattice planes. Adapted with permission from Refs. [174,175]. Copyrights 2012 Institute of Physics and 2011 American Institute of Physics place in the CoSb 3 region, in line with experimental studies. The brittleness of the covalent Sb-Sb bond was compared with the more stable (ionic) Co-Sb bond in interfacial [177] and pure CoSb 3 models [179].
The study of epitaxial SrTiO 3 /SrTi 1−x Nb x O 3 (0 ≤ x ≤ 0.5) modeled by all-electron DFT showed that Nb doping induces metallicity into the originally insulating SrTiO 3 [181]. The Nb doping yields an increase in the carrier concentration while decreasing the thermoelectric response of the heterostructure. The exceptional agreement between the theory and experiment presented in this study demonstrates the high reliability of first-principles calculations for predicting S Coeff within a wide-doping range. In SrTiO 3 interfaces forming epitaxial superlattices with LaNiO 3 (Fig. 11c), DFT + U calculations described how coupling interfaces of opposite polarity generate an electric field in the superlattice [180]. The analyses of the electronic transport and thermoelectric properties showed the anisotropic behavior of this interfacial model, which presents a large positive cross plane S Coeff ~ 135 μV K −1 and a ZT of 0.35 through the p-type superlattice at room temperature.
Understanding the atomic-level interfacial structure, mechanical stability, and electronic features are essential for the design of novel multi-layered materials with an enhanced thermoelectric ZT. These studies showed that the TE efficiency of superlattices could be largely improved compared to their bulk counterparts. The critical parameters for increasing the ZT in superlattices are optimizing the Schottky barrier for enhancing S Coeff while retaining a high σ and reducing the cross plane κ at the interface. The mechanical properties are also an essential factor to consider for engineering applications, since the different thermal expansion coefficient at both sides of the interface may give rise to accumulated thermo-mechanical stresses that detrimentally form interfacial cracks.

Outlook
So far, we have discussed how Boltzmann transport equation, as implemented in BoltzTrap [24,26] can offer solid predictions of the electronic transport properties of the thermoelectric materials. Furthermore, Peierls-Boltzmann transport (PBTE) equation, as implemented in ShengBTE code [136], can determine the phonon lifetimes and κ L . These tools, coupled with the existing material structure databases, either based on computational [183][184][185] or experimental [186,187] compilations, can automate the search for brand new thermoelectric materials through high-throughput calculations within a particular composition or structure space [188,189]. Furthermore, interfaces based on the compounds in these databases can be automatically generated and screened using the MPInterfaces tool [190]. The DFT and transport calculations can be automated using Atomate software package [191] to run automatically, repair, and report the outcome. The outcome then can be searched for suitable intrinsic descriptors such as stability, band gap, flat-band feature, and S Coeff to identify the most promising materials. Recently compounds such as Cu 12 Sb 4 S 13 , Cu 26 V 2 Ge 6 S 32 [192], ZnSnSb 2 , AgInS 2 , AgGaSe 2 , AgInSe 2 , and LiInTe 2 [193] have been predicted to be feasible thermoelectric materials via high-throughput calculations.
High-throughput calculations, although very powerful and predictive, come with some disadvantages. First, the uniformity of the computational settings among many chemical compositions under investigation may cause convergence issues in some compounds. Consequently, the generated data set may not be as accurate for all compounds or even miss some data points for which the DFT calculations catastrophically fail [194]. Moreover, regarding the transport properties, the assumption of a constant relaxation time across a set of compositions may not be realistic and thus is qualitative at best [195]. These problems restrict the capability of most high-throughput calculations in identifying the best thermoelectric materials. Consequently, the promising compounds discovered through such calculations should be further scrutinized with a higher level of theory. For instance, the use of the AMSET package [196] for the electronic transport calculations may remedy some of the shortcomings with the BoltzTrap method as it does not assume a constant relaxation time. AMSET, however, currently relies on a one-dimensional model of the averaged band structure, which limits its predictive power for materials with highly anisotropic band features.
An alternative to high-throughput calculations is to restrict the theoretical search to a handful number of compounds at each investigation-like the studies in Refs [36,39] of Table 1, and Refs. [197,198]-and fine-tune and customize the theoretical settings for the compounds under investigation. As one can see, the balance between the number of compounds investigated and the level of the accuracy is a trade-off to be determined by the constraints on the resources and the design of the search for new thermoelectric materials [199].

Conclusions
This review has highlighted the critical role that advanced computational techniques play in the understanding and design of thermoelectric materials. For instance, band structure calculations based on density functional theory with appropriate functional can aid in predicting the conventional thermoelectric materials operating based on a puddingmold-shaped band near the Fermi level, such as CuAlO 2 , or strongly correlated thermoelectric materials that operate based on Heikes formula for materials with mixed valency such as Na x CoO 2 . Similarly, density functional theory can be employed to study the more complex interface effects at grain boundaries or heterostructures, for instance, in ScN/ ZrN. Furthermore, the wavefunctions calculated by density functional theory can be used to predict the electronic transport (electrical conductivity, Seebeck coefficient, and electronic contribution to thermal conductivity) using Wannier Functions at minimum cost. Phonon dispersion calculations that can elucidate the lattice contribution to the thermal conductivity can be calculated with the Peierls-Boltzmann transport equation. As such, theoretical calculations can predict all properties affecting the thermoelectric performance of a given material or system. Consequently, the computational thermoelectric investigation can complement and assist in the analysis of experimental measurements; for instance, defect concentration, crystal structure, and grain boundary orientations. Moreover, such computational studies can explain the vital atomic-scale features and provide a fundamental understanding of the thermoelectric effect in novel materials for which detailed experimental characterization had not yet completed. Most importantly, such computational investigation can play a predictive role in the design of robust new thermoelectric materials.