Thermoelectric Properties of Sb-S System Compounds from DFT Calculations

By combining density functional theory, quantum theory of atoms in molecules and transport properties calculations, we evaluated the thermoelectric properties of Sb-S system compounds and shed light on their relationships with electronic structures. The results show that, for Sb2S3, the large density of states (DOS) variation induces a large Seebeck coefficient. Taking into account the long-range weak bonds distribution, Sb2S3 should exhibit low lattice thermal conductivity. Therefore, Sb2S3 is promising for thermoelectric applications. The insertion of Be atoms into the Sb2S3 interstitial sites demonstrates the electrical properties and Seebeck coefficient anisotropy and sheds light on the understanding of the role of quasi-one-dimensional structure in the electron transport. The large interstitial sites existing in SbS2 are at the origin of phonons anharmonicity which counteracts the thermal transport. The introduction of Zn and Ga atoms into these interstitial sites could result in an enhancement of all the thermoelectric properties.


Introduction
In a preceding paper [1], we performed a chemical bonding analysis on the ternary Cu-Sb-Se system compounds and showed that the weak interactions, either in local or whole structure, played an important role in lattice thermal conductivity. Since the crystal structure is related to the electronic structure and to the electronic transport path, the analysis of chemical bonds is a bridge between the structure and the thermoelectric properties.
To date, there have been many attempts that tried to explore structure and properties using the concept of chemical bonds, such as the established models of covalent bonding [2], dielectric constants [3], partially ionic binding [4,5] and ionicity [6]. Nonetheless, the complex correlation among physical interactions and bond properties makes the development of a general model very difficult. A commonly accepted concept is that the bonding information can be directly obtained from the charge density. Incontestably, the quantum theory of atoms in molecules (QTAIM) developed by Bader et al. [7][8][9][10] is the most comprehensive density-based topological tool for chemical bonding studies. Bader's quantum theory of atoms in molecules (QTAIM) supplies a well-defined procedure of partitioning a molecule into atomic regions with zero-flux surfaces as boundaries, i.e., the surfaces on which all points satisfy, ∇ρ(r)·n(r) = 0, where n is a unit normal vector of the surface, and then the topology of electron density ρ(r) can be characterized, by analyzing its gradient vector field.
the experimental lattice parameters are not available in the literature. Regarding the SbS 2 compound, we made several structural optimizations on the starting structure extracted from the Materials Project (MP) database (a = b = 6.061 Å, c = 11.59 Å, and V = 426.02 Å 3 ) [30] and found several distinct structures with different lattice constants. We chose to keep the most stable structure (i.e., with the lowest energy, which was about 0.3 eV lower than that from MP) for the subsequent calculations. The Sb 2 S 3 and SbS 2 optimized structures were used as the starting point in the optimization process of the alloyed structures.

Electronic Transport Calculations
Once the electron density was obtained for the optimized structure of interest, the k-point grid was tremendously enlarged to calculate the thermoelectric transport properties. As a rule of thumb, the density was calculated on a grid comprising several thousands (six to ten thousand) k-points. The transport properties were calculated with the BoltzTraP program [35].

QTAIM Calculations
The electron density was analyzed with the program Critic2 [36] that implemented Bader's quantum theory of atoms in molecules (QTAIM). In this work, the model used in Critic2 for the kinetic energy density G(r) calculation was based on the Thomas-Fermi equation with the semiclassical gradient correction proposed by Kirzhnits [37,38], and further refined by the approximate Abramov [39] expression for directly relating G(r) to ρ(r): In turn, the potential energy density was determined using the local Virial equation from the knowledge of G(r) and ∇ 2 ρ(r).
On the basis of the properties of Laplacian distribution and total energy density, the local properties at the BCP could be used to analyze a series of interactions in solids or crystals between pure closed-shell interactions and shared-shell interactions. According to Equation (3), When ∇ρ vanishes at the BCP, G b , V b , and H b are only related to ρ b and ∇ 2 ρ(r) from Equations (4) and (2). Equation (5) displays the variation of H(r)/ρ(r) along with |V b |/G b . The former indicator, the bond degree first introduced by Espinosa [40], stands for the total energy per electron related to the total electronic pressure, and the latter indicator exhibits the competition between potential energy and kinetic energy for bonding formation. These two indicators, which correspond to specific characteristics of bonding, can be used for the bonding classification. At the same time, the kinetic energy per electron G(r)/ρ(r) is related to the variation tendency of H b /ρ b vs. |V b |/G b , which is also an important parameter for classification of bonding [40][41][42].
The crystal structures have been visualized with VESTA [43]. Due to its crystal structure made of one-dimensional (1D) ribbons of polymerized [Sb 4 S 6 ] n , the low-dimensional Sb 2 S 3 is of wide interest in various applications, such as in television cameras [44], microwave devices [45], switching devices, [46] various optoelectronic devices, [47][48][49], as well as photoelectric and thermoelectric cooling devices [50][51][52]. This compound crystallizes into an orthorhombic structure (Pnma) [34], which has two and three different atomic positions for Sb and S, respectively. The Sb1 atoms and Sb2 atoms are located at the centre of a pyramid made by the closest S atoms, a trigonal pyramid for Sb1 and a tetragonal pyramid for Sb2, which together form the infinite [Sb 4 S 6 ] n ribbons parallel to the y-axis (Figure 1a). The [Sb 4 S 6 ] n chains are linked by longer Sb-S interactions and S-S interactions (Figure 1b), forming the whole stable crystal structure. Taking into account the closest and second closest atoms, the coordination numbers are (3 + 4)-fold, (5 + 2)-fold, (3 + 2)-fold, (3 + 2)-fold, and (2 + 6)-fold for Sb1, Sb2, S1, S2, and S3 atoms, respectively. Among the interatomic interactions with the second closest atoms (Figure 2), the Sb1-S1 (3.14 Å), Sb1-S3 (3.20 Å), and S3-S3 (3.58 Å) bonds connect the adjacent ribbons along the z-axis, while the remaining interactions link ribbons in the x direction. Due to these distributions of interatomic interactions, the electronic transport is obviously anisotropic. The electronic motion is preponderant in the y direction and worst in the x one ( Figure S1) but the doping type has a different impact on the electrical conductivity ratio depending on the direction.

Results and Discussion
For Sb 2 S 3 , the experimental energy gap has been reported to lie between 1.61 to 3.8 eV [13,34,[50][51][52][53][54]. The calculated band gaps for bulk Sb 2 S 3 are underestimated relative to these measured values. The value of 1.3 eV, calculated in this work, is in good agreement with the values of 1.29 eV and 1.35 eV determined by GGA-PAW and GGA calculations, respectively [55,56], whereas the local-density approximation (LDA) calculation led to a larger band gap of 1.55 eV [57].
The Engel-Vosko functional and Tran-Blaha modified Becke-Johnson (TB-mBJ) potential can both result in an enlarged band gap that is 1.76 eV for the former and 1.88 eV for the latter [26,58]. In any case, the GGA and LDA calculations are expected to underestimate the gaps [59]. In order to investigate the influence of the band gap on the thermoelectric properties, we performed the thermoelectric properties calculations by using both 1.3 eV and 1.88 eV as band gap energies, the latter one being obtained with the help of a scissor operator. The shapes of the evolution are similar, however, the maximum values of the Seebeck coefficient are different ( Figure 3). In agreement with the results calculated with the TB-mBJ approximation [26], the Seebeck coefficient has maximum values of around ±3000 µV/K near the Fermi level at 300 K when the gap is 1.88 eV, even if, in the present work, we used a denser mesh in the Brillouin zone (3000 vs. 2000 k-points). It can be concluded that the band gap does not influence the thermoelectric properties evolution vs. chemical potential. Only the maximum values of Seebeck coefficient calculated from Boltzmann theory are affected. The underestimation of the band gap results in an underestimation of the Seebeck coefficient that can be corrected by applying a scissor operator. The very high Seebeck coefficient of Sb 2 S 3 originates from the very abrupt density of states (DOS) of valence bands, which are dominated by the Sb-5p and S-3p orbitals [26,[55][56][57][58]. Even if the Sb-S and S-S interactions involving second nearest atoms are weak (see Figure 2), they largely contribute to the DOS near the Fermi level. In effect, the antibonding state of long bonds tends to be filled by electrons and the long range order takes advantage of the largest contributions of S-3p, and also induces astereochemical activity of Sb-5s lone-pair electrons at the top of the valence band. This can be seen from the electronic density of the highest occupied band in Figure 4a. By comparing with Figure 1, we see that S2 and S3 atoms have larger contributions to the electronic density than S1 atoms that have the shortest second nearest neighbor Sb-S bonding, located within the ribbons. The lone-pair electrons on Sb1 atoms have more room than those on Sb2 atoms, which have shorter Sb-S interactions; this indicates that the lone-pair electrons of Sb2 atoms are more stereochemically inert. Moreover, these weaker interactions with delocalized electron state surrounding the ribbons result in the confinement of the dominant electron state within the ribbon and display the quasi-one-dimensional structure of [Sb 4 S 6 ] ( Figure 4). In addition, the soft bonds located close to the (1, 0) point in Figure 2 involve almost all the atoms in the cell, implying a range of low frequency vibrations in the crystal, thus, resulting in low lattice thermal conductivity.
At the same time, n-type doping and p-type doping both display a good thermoelectric performance in the whole temperature range ( Figure S2a), both reaching the optimized thermoelectric power factor around the carrier concentration of ±1 × 10 21 e/cm 3 . However, the n-type thermoelectric performance is severely degraded with further increase in the doping concentration, while the PF of the p-type compound is gradually reduced. The reason is that the Seebeck coefficient has a prominent reduction with the increase in carrier concentration for the n-type compound, whereas the electrical conductivity is substantially enhanced with the increase in carrier concentration for both n-and p-type compounds ( Figure S2b). Irrespective of the carrier concentration up to ±1 × 10 21 cm −3 , it is noticeable that the n-type doping has a larger advantage over the p-type doping. Moreover, the PF exceeds 1 at 300 K and reaches 6.6 at 800 K when doping with 1 × 10 21 e/cm 3 . Taking into account that this compound should have a relatively low thermal conductivity due to soft bonds in the whole cell (see above), Sb 2 S 3 is promising for thermoelectric applications.  Having plotted all the bonds in the cell (Figure 1b), an interstitial site between the ribbons surrounded by a six-membered ring -Sb1-S2-Sb1-S1-Sb2-S3-clearly appears. The packing ratio of the cell is about 29% relative to atomic radii, indicating enough space for introducing impurities. Moreover, the existence of the metastibnite mineral in the Sb-S system, which has the approximate composition of Sb 2 S 3 and contains metals in small amounts [27][28][29], allows one to envisage the adding of impurities into the Sb 2 S 3 crystal. The [Sb 4 S 6 ] n chains, which are parallel to the y-axis, result in a higher electrical conductivity in this direction ( Figure S1). The insertion of metal atoms, in the interstitial sites described above, generates a new chain of atoms along the z-axis that creates infinite [Sb 4 S 6 M 2 ] n chains along this axis by connecting adjacent [Sb 4 S 6 ] ribbons (Figure 5a,b). In order to keep the original [Sb 4 S 6 ] n chains and the semiconducting properties of the material, the light and small Be atom is chosen to be inserted into the cell interstitial sites.

The Be-Sb 2 S 3 Alloy
After structural relaxation, the original Sb 2 S 3 structure is modified, resulting in the formation of pyramidal Be, which connects adjacent [Sb 2 S 7 ] groups and leaves the [Sb 2 S 4 ] groups away from the original [Sb 4 S 6 ] ribbons, as shown in Figure 5c. The alternate layered structures along the x-axis consist of [Sb 2 S 4 ] n chains parallel to the y-axis and infinite [Sb 2 S 9 Be] n chains along the z-axis. According to the topological analysis of chemical bonding (Figure 6), among all the interatomic interactions displayed in Figure 5d, there are four and six different atomic environments for Sb and S atoms, respectively. The Sb1 atoms belong to a square pyramid, the Sb3 and Be atoms belong to a triangular pyramid, and all these polyhedra are contained in the [Sb 2 S 9 Be] group. The [Sb 2 S 4 ] group includes a triangular pyramid containing Sb2 atoms, and Sb4 atoms with two out of the three coordinations, which can be called "bent" coordinations with almost a right angle. In this new cell, in addition to the weaker interactions between secondary neighbors plotted in Figure 6, there are Sb-S and S-S long range interactions (about 5 Å), the BCPs of which are surrounded by the Sb-S nine-membered ring shown in Figure 5d.
The dash line corresponds to the fitting line for Sb-S interactions by least square method.
The energy band structure, given in Figure S3, shows that the two orbitals localized at the top of the valence band almost separate from the other valence orbitals. In addition, the electron density of states of the highest occupied orbital (Figure 7a), as well as that of the top two valence orbitals (Figure 7b), are localized in the [Sb2S4] n chains. As can be seen in Figure 7b, the Sb4-p and S6-p orbitals dominate the electron states of the top two valence orbitals, and the electron density of both Sb4-p and S6-p, which forms a kind of dumbbell, extends along directions so as to form weaker bonds, for example, the Sb4-S6 bonds with a bond length of 3.19 Å, which link the [Sb2S4] n chains along the z-axis. Nevertheless, the electron transport path along the x-axis must pass through the triangular pyramid made of Be, S1, S3, and S5 atoms, which has very poor electron contributions to the DOS near the Fermi level, especially the Be atoms with almost empty states around the Fermi level ( Figure S4). The localized electron states in [Sb2S4] n chains and the electron-poor triangular pyramid both severely restrict the electron transport along the x-axis, resulting in an almost null electrical conductivity in the x direction with p-type doping ( Figure S5a). Since the electron density in the top of the valence band along the yand z-axes is equivalent, without a region with low density, the electrical conductivity is the same in these two directions for p-type doping ( Figure S5a). For n-type doping, the conductivities along yand z-axes are different.
Obviously, the interatomic bonding along the y-axis is stronger than that along the other two axes, where weak bonds between [Sb2S4] n chains and [Sb2S9Be] n chains are located. This is why the electrical conductivity along the y-axis is much higher than that along the other two axes, being similar for electron doping at the bottom of the conduction band.
The Seebeck anisotropy is more complex, and more sensitive to carrier transport and temperature than that of the electrical conductivity; in addition, the scattering mechanisms play an important role [60]. Along the x-axis, the carriers proceed on a path containing weak bonds between layers and pyramidal Be, whereas along the y-axis only the [Sb2S9Be] n chains that stand the pyramidal Be are concerned. Along the z-axis, there are different competing paths, the [Sb2S4] n chains with weak bonds and the [Sb2S9Be] n chains containing pyramidal Be, which should result in different scattering mechanisms for specific carriers and temperature. The Seebeck effect is more prominent along the z-axis and the x-axis for p-type doping and n-type doping, respectively, under the conditions of low carrier concentration and room temperature, as shown in Figure S5b.
Although Sb 2 S 3 Be 2 exhibits an anisotropy in the electronic and thermoelectric properties, the total Seebeck coefficient and the thermoelectric power factor are not as high as those of Sb 2 S 3 ( Figure S6). The main reason is that the electron density at valence band top is confined in the [Sb 2 S 4 ] n chains, making smaller contributions to DOS at the Fermi level ( Figure S7). The extrema of the Seebeck coefficient are −653 µV/K and 595 µV/K at 300 K for the chemical potential of +0.06 eV and −0.065 eV relative to the Fermi level, respectively. These values could be underestimated due to the underestimation of the band gap by GGA calculations. It can be assumed from these extrema values that the n-type doping is more effective than the p-type doping for thermoelectric applications. Calculations of the power factor and the Seebeck coefficient for various carrier concentrations ranging from 5 × 10 18 to 5 × 10 21 cm −3 for both p-and n-type doping ( Figure S8) show, irrespective of the carrier concentration, better properties for electron doping. For carrier concentrations less than or equal to 1 × 10 20 cm −3 , the maximum value of these properties, at a temperature that increases with the carrier concentration, is evidenced. The highest thermoelectric power factor is obtained for the carrier concentration of 1 × 10 21 e/cm 3 , irrespective of the temperature, and its substantial decrease in higher concentrations results from the strong decrease in the Seebeck coefficient. The power factor obtained for 5 × 10 20 e/cm 3 is very close to that for 1 × 10 21 e/cm 3 , for temperatures above 700 K. At lower temperatures, the power factor for 5 × 10 20 e/cm 3 is lower than that for 1 × 10 21 e/cm 3 , this difference being related to higher Seebeck coefficient values below 400 K for 1 × 10 21 e/cm 3 , and especially along the x direction for temperatures below 500 K ( Figure S9).
As for p-type doping, the power factor is lower than that for n-type doping, its highest values being obtained for a carrier concentration of 5 × 10 20 holes/cm 3 in the whole temperature range. The decrease in the power factor for highest carrier concentrations is related to the decrease in the Seebeck coefficient. Moreover, the anisotropy of the Seebeck coefficient for p-type doping is less pronounced than that for n-type doping.

Pure SbS 2
Infinite SbS 2 strings constituted of ψ-SbS 3 tetrahedrons and ψ-SbS 4 bipyramids have been found in BaSb 2 S 4 [61]. SbS 2 has also been synthesized under glass, with a glass transition temperature at 163 • C, and nanocrystalline thin films (Mo:SbS 2 ) forms [62,63]. The crystal structure of SbS 2 has been devised by calculations [30]. SbS 2 has a tetragonal structure and crystallizes in the P4 group. It contains four formula units in the conventional cell for which the optimized lattice parameters are a = 6.56 Å and c = 8.14 Å. This structure exhibits three different positions for the Sb atoms and two for the S atoms, these atoms constituting different [ 4 ], extending to the whole space to form large interstitial sites at 1a and 1d positions, as shown in Figure 8a. A topological analysis of chemical bonds shows that the Sb2 atoms form with S1 atoms, two weak bonds with a bond length of 3.26 Å. The S atoms interact with each other across the interstitial sites, the S2-S2 interactions having longer interatomic distances of 4.33 Å and 5.24 Å. than those of S2-S1 interactions that are equal to 3.52 Å and 3.76 Å (Figure 8b). The electronic band structure ( Figure S10a) calculated with the GGA approach shows that the lowest unoccupied band locates at the A point, whereas the highest occupied band is between the Z and R points, indicating that SbS 2 is an indirect semiconductor with a band gap of 0.66 eV. The relatively low dispersive two top valence orbitals interact only with the other valence orbitals around the G point and are mainly contributed by the Sb2-s and S1-p orbitals and to a lesser extent by the S2-p orbital. The electron density of S1-p and S2-p extends towards the Sb2 atoms and that of Sb2-s towards the S1 atoms (Figure 9a), corresponding to the Sb-S bonds in [Sb2S 4 ] bipyramids and the bond between second nearest atoms Sb2-S1, respectively. The electron states of the two top valence orbitals mainly come from the [Sb2S 4 ] bipyramids, especially from the nearly linear S1-Sb2-S1 fragment bearing an angle of 178.4 • .  The Seebeck coefficient of SbS 2 reaches 1150 and 1025 µV/K at the chemical potential of −0.03 and 0.07 eV relative to the Fermi level at 300 K, respectively ( Figure S11a), indicating that a p-type doping is more efficient for TE applications. These results are in agreement with a larger DOS derivative observed near the Fermi level in the valence band, which is contributed by the [Sb2S 4 ] bipyramids. The single peak of the electrical conductivity observed at positive chemical potential near the Fermi level and related to the single peak observed in the conduction part of the DOS, which originates from the [Sb1S 4 ] tetrahedron, weakens the thermoelectric effect. As a consequence, the thermoelectric power factor is poorer for n-type doping than for p-type doping. Even at 300 K in the chemical potential range [−0.03 Ry, 0.04 Ry], which is larger than the usual doping range of~10 k B ·T, the conductivity is the same along the x, y, and z directions ( Figure S11b), an anisotropy is observed for the electrical conductivity in the remaining chemical potential ranges, in concordance with the crystal structure. This anisotropy can be explained by the electron transport path that is different along the x, y and z directions. The electrons must go through the six-membered polyhedral rings along the z-axis, whereas they have two optional paths along the xand the y-axes. Indeed, the electrons can go through the [Sb1S 4 ] tetrahedron and the [Sb2S 4 ] bipyramid or the [Sb3S 4 ] tetrahedron and the [Sb2S 4 ] bipyramid. Obviously, the more efficient electron path in the cell is that passing through the [Sb1S 4 ] tetrahedron and the [Sb2S 4 ] bipyramid, corresponding to better electrical properties along the xand the y-axes than along the z-axis.
To be more accurate, the thermoelectric properties investigation has to be completed by considering the influence of the doping. Although the Seebeck coefficient has similar evolution vs. temperature irrespective of the doping type, the p-type doping has a larger Seebeck coefficient than that of the n-type doping for the same carrier concentration ( Figure S12a).
The Seebeck coefficient decreases with the increase in carrier concentration and behaves distinctively with respect to temperature depending on the carrier concentration. For low carrier concentrations, it increases up to around 500 K, and then drastically decreases above this temperature, whereas, for high carrier concentrations, it slightly increases with temperature to reach a plateau at high temperatures. Regarding the thermoelectric power factor (Figure S12b), the influence of the temperature differs from n-type to p-type doping.
The p-type doping yields a higher thermoelectric power factor at low temperatures, whereas the n-type doping performs better at high temperatures, the crossing point being displaced towards higher temperatures as the carrier concentration increases. For the hole concentration of 1 × 10 20 cm −3 , the power factor reaches its maximum at about 400 K, whereas for the same electron concentration, the power factor reaches its maximum at about 850 K. For higher carrier concentrations, the p-type doping shows better efficiency than the n-type doping in the whole temperature range of interest (0-1000 K). For both p-type and n-type doping, the maximum power factor is obtained for a carrier concentration of 5 × 10 20 cm −3 . Although the thermoelectric performance of SbS 2 seems to be not as good as that of Sb 2 S 3 , the large interstitial site in the cell and the long-range S-S interactions should be at the origin of phonon anharmonicity and should result in low lattice thermal conductivity.
The packing ratio of the SbS 2 cell is about 27% relative to atomic radii, indicating that interstitial sites at 1a and 1d Wyckoff positions are large enough to accommodate foreign metal atoms.
Two situations can occur when introducing metal atoms in these interstitial sites. The first situation occurs when the metal atoms form strong bonds with adjacent S atoms inducing a decrease in the cell volume. The second situation occurs when the metal atoms form weak interactions with adjacent S atoms and result in an increase in cell volume. The latter situation is similar to a cage structure, making the metal atoms vibrate in the cage, which greatly enhances the vibrations anharmonicity, and thus reduces the lattice thermal conductivity. Both situations are discussed hereafter.

The Zn-SbS 2 Alloy
The first situation can be achieved by introducing two Zn atoms at 1a and 1d positions. After structure optimization, the cell volume of SbS 2 Zn 2 decreases by about 6.7%, resulting from an increase in a and a decrease in c parameters.
The energy densities of all the Sb-S bonds in SbS 2 and SbS 2 Zn 2 are plotted in Figure 10. The bond lengths in the [Sb1S 4 ] and [Sb3S 4 ] tetrahedrons increase as the coordinated atoms of Sb1 and Sb3 change from S2 to S1 and from S1 to S2, respectively, when the composition goes from SbS 2 to SbS 2 Zn 2 . At the same time, the bond angle defined by S1-Sb2-S1 decreases from 178.4 • to 154 • and the corresponding bond lengths decrease.
By contrast, the bond angle defined by S2-Sb2-S2 increases while the corresponding Sb2-S2 bond lengths remain the same, making the [Sb2S 4 ] bipyramids more planar, as shown in Figure 11a. In spite of the formation of strong Zn-S bonds in the [ZnS 4 ] tetrahedrons with a bond length of 2.36 Å, which is slightly longer than the sum of the atomic radii (2.3 Å), the long-range S1-S2 interactions still exist, although with elongated interactomic distances (5.69 Å), as shown in Figure 11b. With two Zn atoms introduced at 1a and 1d Wyckoff positions, under GGA calculations, the band gap decreases from 0.66 to 0.52 eV. Some very localized bands appear in the energy range of −6 to −7 eV relative to the Fermi level ( Figure S13a). These bands originate from the Zn-3d orbitals and generate a very intense DOS ( Figure S13b and inset), indicating that Zn-S bonds are strong and tend to lower energy. The maximum of the highest occupied orbital, which is dominated by Sb1-s and S1-p (Figure 12a), locates at the G point. The two top valence orbitals almost separate from other bands along the high symmetry k-points' path ( Figure S13a), the electron states of which are dominated by the S1-p and S2-p orbitals (Figure 12b). The hole conduction properties are governed by these two orbitals ( Figure S14a). Irrespective of the doping type, the electronic transport properties along the xand y-axes are the same, whereas they differ from those along the z direction, except in the chemical potential range 0-0.68 eV. Along the z-axis, the optimal conduction path goes through the [ZnS1 4 ] and [Sb1S1 4 ] polyhedra, whereas along the xand y-axes, it goes through a path alternating S1 and S2 atoms. Thus, the highest electrical conductivity observed along the z-axis for hole doping can be attributed to the dominance of S2 atoms at the top of the valence band. Regarding the thermoelectric effect, the very low density of states at the valence band edge ( Figure S13b) results in a Seebeck coefficient for hole doping lower than that for electron doping ( Figure S14b). Nonetheless, the Seebeck coefficient reaches 797 and 920 µV/K at the chemical potential of −0.06 and 0.04 eV, respectively. The higher Seebeck coefficient for n-type doping is in agreement with the thermoelectric power factor calculated for various doping levels ( Figure S15a). The maximum thermoelectric power factor for hole doping with the optimal carrier concentration of 1 × 10 21 h/cm 3 is even much lower than that with electron doping with a carrier concentration of 1 × 10 20 e/cm 3 . However, the Seebeck coefficient for electron doping, at a low carrier concentration, drastically falls at high temperatures, as seen in Figure S15b. Overall, SbS 2 has a good thermoelectric power factor for electron doping when introducing Zn atoms in the interstitial sites but the anharmonicity sources present in SbS 2 , namely the long range S1-S2 interactions and the secondary Sb-S interactions, are weakened.

The Ga-SbS 2 Alloy
When two Ga atoms are introduced into the 1a and 1d interstitial sites of SbS 2 , the cage structure after optimization is maintained. The volume of the new cell is increased by about 10% as the a lattice parameter and the c lattice parameter increase and decrease, respectively. The Ga atoms are surrounded by eight adjacent S atoms forming weak Ga-S interactions with interatomic distances from 3.14 to 3.38 Å (Figure 13), which are enclosed in the range delineated by the atomic radii sum (2.24 Å) and the Van der Waals radii sum (3.67 Å). Meanwhile, the interatomic distance for the weak S1-S2 interactions still reaches 3.6 Å, which is close to the sum of the Van der Waals radii. The architecture around the Sb atoms has undergone only slight changes, except that the bond length of Sb1-S2 increases by about 0.22 Å. In the band structure, the top valence orbital separates from other valence orbitals and, as observed for SbS 2 , the bottom conduction orbital is completely independent from other conduction orbitals ( Figure S16). Figures 14 and 15 show that the electron states of the top valence orbital and the bottom conduction orbital are dominated by Sb1-s, S2-p and Sb3-s, S1-p, respectively. These electron states belong to the [Sb1S 4 ] and [Sb3S 4 ] tetrahedrons, respectively, which constitute the specific electron transport path from the valence band to the conduction band.   The Seebeck coefficient for n-type doping is larger than that for p-type doping, as shown in Figure S17a. It reaches 525 and −853 µV/K at 300 K, for the chemical potential of −0.1 and 0.001 eV, respectively. Regarding the thermoelectric power factor (Figure S17b), for similar carrier concentrations up to 10 20 cm −3 , at low temperatures (<500 K), the PF is higher for n-type doping. For higher temperatures, the PF is higher for p-type doping irrespective of the concentrations. For both p-type and n-type doping, the PF increases with the carrier concentration, up to 5 × 10 20 cm −3 and decreases beyond this value. As with SbS 2 , along the z-axis, the electrons must go through the six-membered polyhedral rings, whereas along both the xand y-axes they have two possible paths depending on the doping type, namely through the [Sb1S 4 ] tetrahedrons and [Sb2S 4 ] bipyramids in the valence band or through the [Sb3S 4 ] tetrahedrons and the [Sb2S 4 ] bipyramids in the conduction band. Thus, the electrical properties along the xand y-axes are similar and better than along the z-axis ( Figure  S18a). The anisotropy between the x/y and z directions is less pronounced for the Seebeck coefficient, except for the highest energy levels in the valence band ( Figure S18b).

Concluding Summary
By combining the analysis of the distribution of chemical bonds through the QTAIM approach, with that of the electronic band structure, density of states, and electron density maps, a detailed description of the structure-properties relationships can be drawn. The electronic transport paths going through specific atom groups have been distinctly identified in the cell, being responsible for the isotropic or anisotropic character of the electrical properties. As for the Seebeck coefficient, the isotropy or anisotropy are much more difficult to explain, as the Seebeck is comprehensively influenced by carriers, temperature, scattering mechanisms, and lattice, which need to be discussed with specific doping.
In Sb 2 S 3 , the electronic transport is anisotropic. The Sb electron lone pairs do not have the same activity, being more or less stereochemically active, depending on their environment. Together with the fact that soft bonds involve almost all the atoms in the cell, the lattice thermal conductivity should be low.
In Sb 2 S 3 Be 2 , the structure and density of states analysis shows that, the electron transport in the x direction has to pass through the Be atom arranged in triangular pyramid with the neighboring S atoms, both showing low electronic density near the Fermi level, thus, resulting in nearly null electrical conductivity. This situation seems to be more favorable for the Seebeck effect.
The specific arrangement in bipyramids and tetrahedrons of the antimony atoms in SbS 2 favors the electron transport in the x and y direction of the cell. Alloying SbS 2 with Zn results in a 6.7% decrease in cell volume. Along the z direction, the good electron conduction is attributed to the presence of ZnS 4 and SbS 4 polyhedra and to the dominance of the p states of sulphur atoms at the top of the valence band. Overall, the electronic thermoelectric properties of this alloy are good, but the lattice thermal conductivity is predicted to be higher than that in pure SbS 2 . As SbS 2 is alloyed with gallium, the cell volume increases by about 10%. The arrangement around the antimony atoms is weakly affected. A quite specific band structure is exhibited by this structure with both a single valence band and a single conduction band around the Fermi level, well separated from the other valence and conduction bands. These isolated bands are both contributed by Sb-s and S-p orbitals and the SbS 4 atomic arrangements constitute the corresponding transport path for electrons from the valence band to the conduction band. As it turns out, the power factor is lower for SbS 2 alloyed with gallium than for SbS 2 alloyed with zinc. All in all, the Sb-S compounds show intermediate thermoelectric performances that can be compensated for by their environmental friendliness.
Supplementary Materials: The following are available online at http://www.mdpi.com/1996-1944/13/21/4707/s1, Figure S1: Electrical conductivity as a function of chemical potential in different directions at 300 K for Sb2S3. Figure S2: Thermoelectric properties vs. temperature for Sb2S3 with various carrier concentrations (e/cm3 and holes.cm-3). a) Power factor; b) Electrical conductivity. Negative and positive values of carrier concentrations correspond to n-type and p-type doping, respectively. Figure S3: Energy band structure of Sb2S3Be2. Figure S4: Density of states projected on the atoms involved in the Be pyramid: a) Be atom; b,c,d) S atoms. Figure S5: (a) Electrical conductivity; (b) Seebeck coefficient in different directions as a function of chemical potential for Sb2S3Be2 at 300 K. Figure S6: Thermoelectric properties (Seebeck coefficient (a), electrical conductivity (b) and power factor (c)) as a function of the chemical potential at 300 K for Sb2S3Be2. Figure S7: Density of states (DOS) of Sb2S3Be2. Figure S8: Thermoelectric properties of Sb2S3Be2: (a) Thermoelectric power factor; (b) Seebeck coefficient as a function of temperature for various carrier concentrations (cm-3). Figure S9: Seebeck coefficient of Sb2S3Be2 vs. temperature for various carrier concentrations (cm-3) and directions. Figure S10: (a) Electronic band structure of SbS2; (b) Density of states (DOS) of SbS2. Figure S11: (a) Thermoelectric properties of SbS2 and (b) Electrical conductivity along the different directions (x,y,z) as a function of chemical potential at 300 K. Figure S12: Thermoelectric properties of SbS2 vs. temperature for various carrier concentrations (cm-3) (a) Seebeck coefficient and (b) thermoelectric power factor. Figure S13: (a) Band structure of SbS2Zn2 and (b) Density of states (DOS) of SbS2Zn2.The insert in (b) corresponds to the PDOS of Zn atoms. Figure S14. Electronic properties of SbS2Zn2 (a) electrical conductivity for x,y and z directions; (b) thermoelectric properties vs. chemical potential at 300 K. Figure S15: Thermoelectric properties of SbS2Zn2: (a) Thermoelectric power factor and (b) Seebeck coefficient vs. temperature for various carrier concentrations (e/cm3). Figure S16: (a) Band structure of SbS2Ga2 and (b) Density of states (DOS) of SbS2Ga2. Figure S17: Electronic properties of SbS2Ga2 (a) Thermoelectric properties vs. chemical potential at 300 K; (b) Power factor as a function of temperature for various carrier concentrations (cm-3). Figure