Strain Effects on the Electronic and Thermoelectric Properties of n(PbTe)-m(Bi2Te3) System Compounds

Owing to their low lattice thermal conductivity, many compounds of the n(PbTe)-m(Bi2Te3) homologous series have been reported in the literature with thermoelectric (TE) properties that still need improvement. For this purpose, in this work, we have implemented the band engineering approach by applying biaxial tensile and compressive strains using the density functional theory (DFT) on various compounds of this series, namely Bi2Te3, PbBi2Te4, PbBi4Te7 and Pb2Bi2Te5. All the fully relaxed Bi2Te3, PbBi2Te4, PbBi4Te7 and Pb2Bi2Te5 compounds are narrow band-gap semiconductors. When applying strains, a semiconductor-to-metal transition occurs for all the compounds. Within the range of open-gap, the electrical conductivity decreases as the compressive strain increases. We also found that compressive strains cause larger Seebeck coefficients than tensile ones, with the maximum Seebeck coefficient being located at −2%, −6%, −3% and 0% strain for p-type Bi2Te3, PbBi2Te4, PbBi4Te7 and Pb2Bi2Te5, respectively. The use of the quantum theory of atoms in molecules (QTAIM) as a complementary tool has shown that the van der Waals interactions located between the structure slabs evolve with strains as well as the topological properties of Bi2Te3 and PbBi2Te4. This study shows that the TE performance of the n(PbTe)-m(Bi2Te3) compounds is modified under strains.


Introduction
Energy production and environment protection have become two of the most critical current issues, as fossil fuels constitute a finite source of energy, and at the same time their consumption drastically affects the climate patterns. Only a fraction of the energy released by the burning of fossil fuels is converted into mechanical energy or electricity, whereas most of that energy is released as heat. Therefore, policy makers and scientists alike are exploring new routes to provide green, clean and efficient energy. Thermoelectric generators (TEGs) are solid-state devices that convert a heat flux directly into electrical power [1]. The heat-to-electricity conversion exploits the Seebeck effect. A thermoelectric (TE) converter is constituted of n-type and p-type materials, and the absence of moving or mechanical parts allows them to function for a substantially long time without the need for repair. In addition, they are extremely quiet, reliable and scalable, making them ideal for small distributing power generation. One of the main drawbacks of the TE converters is the low conversion efficiency of their constitutive TE materials [2].
In the past decades, many TE materials systems have been studied, such as Half-Heusler alloys [3][4][5][6], chalcogenides [7][8][9][10], skutterudites [11][12][13][14], clathrates [15][16][17][18] and zintl phases [19][20][21][22]. The performance of TE materials is determined by the dimensionless figure of merit zT = S 2 σ/(κ e + κ l ), where S is the Seebeck coefficient, σ is the electrical conductivity and κ e and κ l are the electronic part and lattice part of thermal conductivity. High TE efficiency needs high S and σ and low thermal conductivity, κ. Improving TE performance has been a big challenge. To meet this challenge, various approaches that can be classified into two categories, namely band engineering and phonon engineering, have been used. Through constructing superlattice structures or creating nanostructures, the phonon-boundary scattering can be enhanced, resulting in an extremely low thermal conductivity [23,24]. On the other hand, the TE properties can be boosted via enhancing effective density of states (DOS) by band convergence near the Fermi level [2,25,26] or changing the forbidden bandwidth by strains [26]. In the present work, tensile and compressive biaxial strains have been applied to Bi 2 Te 3 , PbBi 2 Te 4 , PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 to investigate their effects on the thermoelectric and bonding properties of these compounds. The last section is devoted to the investigation of the topological insulator properties of Bi 2 Te 3 and PbBi 2 Te 4 under peculiar strains.

Methods and Computational Details
The calculations of structural, electronic and thermoelectric properties have been performed within the frame of density functional theory using the all-electron FP-LAPW approach with the local orbital method, as implemented in WIEN2K (version 19.1, 2019, Technology University of Vienna, Vienna, Austria) [27]. Several exchange-correlation functionals have been used, the details of which will be mentioned when appropriate in the results section. Although hybrid exchange-correlation functionals usually yield excellent results compared with experimental ones regarding band-gap energies, as exemplified in [28,29], the large number of calculations performed in the present work preclude the use of such time-consuming functionals. Therefore, the results presented in this work have been obtained with pure, density-based functionals. For structural optimizations, the Brillouin zone has been sampled with the k-meshes 8 × 8 × 8, 8 × 8 × 8, 12 × 12 × 2 and 12 × 12 × 2, for Bi 2 Te 3 , PbBi 2 Te 4 , PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 , respectively. For the subsequent convergence of the self-consistent energy, the Monkhorst-Pack k-meshes have been set as 16 × 16 × 16, 16 × 16 × 16, 18 × 18 × 4 and 18 × 18 × 2. The valence electrons for Pb, Bi and Te have been taken as 5d 10 6s 2 6p 2 , 5d 10 6s 2 6p 3 and 4d 10 5s 2 5p 6 . The total energy and atomic forces' convergence thresholds have been defined as 0.136 meV and 0.257 meV/Å, respectively. The RmtKmax value has been set to 9.0 and the Radius of Muffin Tin (RMT) used for both Bi and Te atoms in this study has been set to 2.5 Å. The structure and charge density calculated above are used to analyze the topological properties within QTAIM theory [30]. The transport properties of the compounds with and without strain have been calculated with the BoltzTraP2 code [31] based on the use of a full band's structure in the Brillouin zone. The sampling, which is important in transport calculation, has been performed with a very dense k-mesh of 36  The surface states of Bi 2 Te 3 and PbBi 2 Te 4 , have been calculated from 2D film structure made of six quintuple-layer slabs and six septuple-layer slabs, respectively. The optimization of the films has been performed with a vacuum height of 15 Å to avoid the artificial interaction between atom layers. The k-meshes used to sample the Brillouin zone have been set to 10 × 10 × 1. The parameters of the total energy and atomic forces' convergence, RmtKmax and RMT values are identical to those used in the self-consistent calculations of the bulk compounds.

Compounds' Structural Information
Bulk Bi 2 Te 3 and PbBi 2 Te 4 crystallize in the rhombohedral lattice system (R3m space group) with 5 and 7 atoms in the primitive cell stacked along the c-axis, respectively. However, they can also be described with a hexagonal cell constituted by three 5-atomlayered slabs and three 7-atom-layered slabs for Bi 2 Te 3 and PbBi 2 Te 4 , respectively ( Figure 1). The PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 compounds have a hexagonal unit cell and belong to the P3m1 space group. The crystal structure of PbBi 4 Te 7 is a 12-atom-layered one consisting of a 5-atom-layered slab and a 7-atom-layered one. Likewise, Pb 2 Bi 2 Te 5 is constituted by one 9-atom-layered slab. Two possible atom sequences were found for Pb 2 Bi 2 Te 5 , one by Petrov [32] (-Te-Pb-Te-Bi-Te-Bi-Te-Pb-Te-) and the other one by Chatterjee [33] (-Te-Bi-Te-Pb-Te-Pb-Te-Bi-Te-). In the present study, we used the most stable sequence as found by Ma et al. [34], which corresponds to the Chatterjee one. Pb 2 Bi 2 Te 5 with this sequence was found to be a semi-conducting compound. In all of these structures, the slabs are linked together by Te-Te interactions. The exchange-correlation functionals we used for structure optimization were the Local Density Approximation (LDA) [35,36], the Perdew-Burke-Ernzerhof (PBE) [37] as well as the rev-vdw-DF2 [38] to account for Van der Waals forces that have to be considered in these compounds due to the weak Te-Te interactions between the slabs of the n(PbTe)m(Bi 2 Te 3 ) system. The optimized lattice parameters have been determined from total energy minimization with respect to the crystal cell volume and c/a ratio. The fitting result was obtained by using the Birch-Murnaghan equation of state: where V 0 and V represent the initial and deformed volume respectively, B 0 is the bulk modulus and B P is the derivative of the bulk modulus with respect to pressure. Based on the optimized lattice constants, the total residual force on all atoms was relaxed to less than the aforementioned threshold. The values of the equilibrium lattice parameters and atomic positions in each of the compounds are provided in Table 1 together with experimental data from the literature. Starting from the optimized structure, strains in the range −10% to 10%, which could experimentally originate from epitaxial strains in thin films [42], have been applied to a and b (b = a) lattice parameters. Based on the calculated mechanical properties [34] of the Pb 2 Bi 2 Te 5 compound, we assume that all the compounds of interest can withstand such strains. Indeed, the bulk modulus, shear modulus and Young modulus amount to 45.6, 32.0 and 77.9 GPa respectively, for Pb 2 Bi 2 Te 5 . The strain η is calculated as η = (a − a 0 )/a 0 . Then, the interlayer distances have been optimized for all 4 compounds under the applied strain. Figure 2 shows the calculated inter-slab distance (denoted slab distance) and the Te-Te distance between adjacent slabs with respect to strain (compressive strain when η < 0, tensile strain when η > 0). Under an increasing compressive strain, slab distance and Te-Te distance between adjacent slabs increase simultaneously for all four compounds, with the latter increasing more rapidly than the former. Further, for PbBi 2 Te 4 and Pb 2 Bi 2 Te 5 , both the Te-Te distance and the slab distance increase more rapidly than for Bi 2 Te 3 and PbBi 4 Te 7 , indicating lower van der Waals interactions between the Te atoms. By contrast, under an increasing tensile strain, slab distance and Te-Te distance between adjacent slabs decrease simultaneously for all four compounds, but for all the compounds, the Te-Te distance decreases less than the slab distance.   Figure 3b). The reference energy is taken as the energy of the structures under 10% compressive strain, which corresponds in all cases to the highest value. The SOC effect has been investigated owing to the presence of Pb and Bi heavy metal elements in the n(PbTe)m(Bi 2 Te 3 ) compounds. One can see that for all compounds, with and without SOC, the minimum energy corresponds to a very low strain (η < 2%), testifying to the stability of the unstrained structure. The SOC effect is very weak, and the biggest difference is observed for the relative energy of Bi 2 Te 3 and PbBi 2 Te 4 under tensile strains. By contrast, as shown in our previous study [34], the SOC effects are non-negligible on both the electronic and the thermoelectric properties. Therefore, the SOC effect has been accounted for in the forthcoming calculations.

Electronic Structure
The calculated DOS of Bi 2 Te 3 , PbBi 2 Te 4 , PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 with and without strain are depicted in Figure 4. On the basis of the conduction band minimum (CBM) and valence band maximum (VBM) positions, all the strain-free compounds are found to be indirect semiconductors. The VBM of Bi 2 Te 3 and PbBi 2 Te 4 locates along the L-Z path, while the CBM locates along the Z-F path (see Supplementary Figures S1 and S2), which is in agreement with previously reported data [43]. Under 3% compressive strain, irrespective of the compounds, the DOS show a semiconducting behavior. At 6%, all the compounds are semiconductors except Bi 2 Te 3 , and at 9%, they are all metals except PbBi 2 Te 4 . Under tensile strains, all the compounds are metals except PbBi 2 Te 4 at 3%. To give a clear evolution of the semiconductor-to-metal transition, we have precisely calculated the CBM and VBM energy values from the DOS using a very dense k-point mesh. The CBM and VBM energy evolution as a function of strain is depicted in Figure 5, as well as the energy gap (E CBM − E VBM ).
In absence of strain, Bi 2 Te 3 , PbBi 2 Te 4 , PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 are narrow semiconductors. The compounds remain semiconductors in a wide range of strain: from −4% to 0% for Bi 2 Te 3 , from −9% to −3% and −1% to +3% for PbBi 2 Te 4 , from −7% to 0% for PbBi 4 Te 7 and −6% to +1% for Pb 2 Bi 2 Te 5 . As it turns out, the range of non-zero band gaps is larger under compressive strains than under tensile ones. Insets in Figure 5a, b show the evolution of VBM and CBM energy versus strains at Z point for PbBi 2 Te 4 and at Γ point for Bi 2 Te 3 . The energy gap of PbBi 2 Te 4 with 3% compressive strain is very small and direct at Z point, whereas it is indirect for most of the strained compounds. The energy difference at the Γ point for Bi 2 Te 3 vanishes for a compressive strain of about 8%.

Transport Properties
Based on the calculated electronic structure, the Seebeck coefficient (S) of Bi 2 Te 3 , PbBi 2 Te 4 , PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 has been investigated. Figure 6 shows the calculated S at 300 K as a function of doping level under different strains for all the compounds. We found that the maximum Seebeck coefficients for p-type Bi 2 Te 3 , PbBi 2 Te 4 , PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 are located at −2%, −6%, −3% and 0% strain, which are near the maximum band gaps of each compound. For Bi 2 Te 3 , with increasing compressive strain up to −2%, the Seebeck coefficient increases for both p-type and n-type doping. For higher compressive strains, the Seebeck coefficient decreases while remaining higher than that of the strain-free compound up to −3% (p-type) and −4% (n-type). The maximum Seebeck coefficient for PbBi 2 Te 4 is obtained at a compressive strain of −6% for p-type doping and −8% for n-type doping. In order to explain the difference in the Seebeck coefficients under strain, the tight relationship between S and electronic structure has been studied. Following the band theory, the Seebeck coefficient can be described as [44,45]: where k B is the Boltzmann constant, e is the elementary charge, E F and E G are the Fermi energy and the band-gap energy, N v and N c are the effective DOS of the valence and conduction bands, µ n and µ p are the mobility of electrons and holes and n and p are the number of electrons and holes, respectively.
The Seebeck coefficient depends on the carrier concentration, the effective DOS of the conduction and valence band, band-gap width and Fermi energy. By means of strainmodified DOS, the Seebeck coefficient can be changed accordingly. The maximum Seebeck coefficients at 300 K of PbBi 2 Te 4 without strain are 198 µV/K and −81 µV/K for p-type and n-type doping, respectively. These results are in agreement with the effective DOS values and the band gap, i.e., the value of the valence band is higher than that of the conduction one (Figure 4b). Besides, due to the similar effective DOS of the valence band and conduction band of Bi 2 Te 3 , the Seebeck coefficient magnitudes for p-type and n-type doping are very close. As can be seen in Figure 4d and Supplementary Figure S2  For solving the Boltzmann transport equation, the BoltzTraP2 code uses the relaxation time approximation, and therefore the scattering time τ remains undefined. In general, for a given structure, τ depends on temperature and doping level. We can obtain a temperature-dependent τ from the so-called deformation potential theory based on the deformation potential, the carriers' effective mass and the elastic constants. Nevertheless, with the changes of strain, the effective mass of carriers also changes, which would make it computationally difficult to evaluate the scattering time for each and every case. Therefore, we have kept τ as is in the electronic conductivity σ. Due to the low Seebeck coefficient in metallic compounds, only the electronic conductivity results for strained compounds with opened band gap are reported in the following. The corresponding strain ranges are −5% to 0% for Bi 2 Te 3 and PbBi 4 Te 7 , −8% to 2% for PbBi 2 Te 4 and −4% to 1% for Pb 2 Bi 2 Te 5 . Figure 7 shows the calculated σ/τ at 300 K as a function of doping level for various strains. In most cases, the electrical conductivity decreases with the increase of compressive strain for both p-type and n-type doping, at least for high doping levels. The peculiar behavior of Pb 2 Bi 2 Te 5 can be noticed for n-type doping as the electrical conductivity is rather insensitive to the applied strains. Except for Bi 2 Te 3 , which shows similar evolutions of the electrical conductivity under strains for both pand n-type doping, the electrical conductivity is more sensitive to the applied strain for p-type doping than for n-type doping. This is probably due to the larger dispersion of the conduction bands of PbBi 2 Te 4 , PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 than that of the valence bands around the Fermi level (see Figure 4) in the range of strains with opened gaps. According to the transport equations, S and σ evolve in opposite ways from one another. This is still valid for our results under strains, which indicates that strains could not de-correlate these quantities. Hence, the power factor, PF = S 2 σ, is essential to evidence the evolution of the material efficiency with respect to the applied strains. The PF of the investigated p-type and n-type compounds is depicted in Figure 8 for various compressive and tensile strains that have been chosen in order to better evidence the evolution of the power factor. Overall, the power factor is improved by applying strains. Specifically, the power factor for p-type compounds tends to increase upon tensile strains, whereas that of the n-type compounds tends to increase under compressive strains. The largest PF values are obtained for the p-type compounds under tensile strains.
The calculated power factors are fairly good if one considers τ to amount to a commonly used value of 10 −14 s [46][47][48]. Some of these values are, however, calculated for closed energy gaps, which highlights the difficulty to optimize the thermoelectric properties of these materials.

Charge and Bonding Evolution under Strains
In order to complete the investigation, we have calculated charge density differences between strained and unstrained structures (ρ diff = ρ unstrained − ρ strained ). Table 1 shows a great difference in lattice parameters and atomic positions between the strained and the unstrained structures. In order to compare the difference of charge density between strained and unstrained structures, we have meshed the cells in the [001] direction with grids leading to the same number of ac-planes (00x) between each pair of adjacent atoms. Figure 9 displays the charge differences in Bi 2 Te 3 , PbBi 2 Te 4 , PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 for compressive (−10% and −5%) and tensile strains (5% and 10%) with an iso-surface of 0.05 e/bohr 3 . The yellow color represents charge increase, while the cyan one corresponds to charge decrease. A compressive strain leads to charge increase in the in-plane (ab-axis) direction, and decrease in the cross-plane (c-axis) direction, with the decrease being larger for a higher value of the compressive strain (cf. arrow in Figure 9). Conversely, opposite results are obtained under tensile strain, which reveals that the charge density between the slabs increases when going from a compressive strain to a tensile one. As a consequence, the van der Waals bonds located between the slabs should be strengthened under a tensile strain.
To further understand the evolution of the bonds under strains, we have used the quantum theory of atoms in molecules (QTAIM) [49]. All irreducible bond critical points (BCP) in Bi 2 Te 3 and PbBi 2 Te 4 are considered. Due to similar properties between b2 and b5, and b3 and b6 in PbBi 4 Te 7 (see Supplementary Figure S5), and b4 and b5 in Pb 2 Bi 2 Te 5 (see Supplementary Figure S6), only b1 to b4 BCPs are taken into account for PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 . The charge density, ρ, and its Laplacian, ∇ 2 ρ, at these BCPs, presented in Figure 10, reveal that the strain affects the bonds between slabs much more than the in-slab ones. As shown by the positive values of the electron density Laplacian irrespective of the strain, there is no charge accumulation at the considered BCPs. The values of the Laplacian, ∇ 2 ρ, also called valence shell charge concentrations (VSCC), and the charge density, ρ, at the in-slab BCPs diminish when going from a compressive strain to a tensile one, indicating that the in-slab bonds become weaker, while the interaction between slabs is getting stronger. These results are in agreement with those related to the charge density differences presented above.  Based on Bader's theory of Atoms in Molecules, the sign of ∇ 2 ρ can be used to differentiate the closed-shell (CS) interactions from the shared-shell (SS) interactions. Closed-shell interactions (ionic, H-bonds and vdW) have a large positive value of ∇ 2 ρ, |V|/G < 1 and a small ρ. Conversely, ∇ 2 ρ < 0, |V|/G > 2 and a large ρ are expected for shared interactions (covalent or polar bonds). In our study, all BCPs, except the Te-Te ones, show a small value of ρ, positive ∇ 2 ρ with 1 < |V|/G < 2 ( Figure 11). As a consequence, the corresponding bonds cannot be considered as pure covalent nor pure closed-shell ones. They belong to a transit region identified by Espinosa [50] and confirmed by Dinda [51]. By contrast, as expected, the Te-Te bonds with a small value of ρ, positive ∇ 2 ρ with |V|/G < 1 correspond to pure closed-shell interactions. When going from compressive to tensile strain, the bond degree defined by Espinosa as BD = H/ρ goes from a positive value to a negative one, leading, as already mentioned above, to a bond strengthening for the van der Waals bonds located between the slabs.  [52][53][54]. As discussed above, band structure changes along with the strain. In addition, several Dirac-like cones have been observed along the high symmetry path. Based on the results shown in Figure 5a,b, two interesting strained electronic structures, presenting gapless bands and Dirac-like dispersion, have been selected to investigate the topological properties of Bi 2 Te 3 and PbBi 2 Te 4 . The corresponding strain values are −8.4% and −2.2% for Bi 2 Te 3 and PbBi 2 Te 4 , respectively. All the results presented hereafter have been obtained with these strains.
The band structures and partial density of states (PDOS) are depicted in Figure 12. The electronic gaps are closed at the Γ and Z points for Bi 2 Te 3 and PbBi 2 Te 4 , respectively. The PDOS indicates that the main contributing orbitals close to the Fermi level are Te 5p, Bi 6p and Bi 6s for both Bi 2 Te 3 and PbBi 2 Te 4 . For PbBi 2 Te 4 , Pb 6s orbitals also significantly contribute. More specifically, the valence band of Bi 2 Te 3 near Γ is composed predominantly of Te states, whereas the conduction band is contributed by Bi states. A similar situation is found at Z point for PbBi 2 Te 4 . Such an inversion of gap edges may indicate a change of the parity of the occupied states, revealing a promising topological insulator (TI) [55]. In order to investigate the surfaces' electronic states, thin film structures of Bi 2 Te 3 and PbBi 2 Te 4 , both with a P3m1 symmetry, have been designed. These structures, which consist of six 5-atom-layered slabs (5L-5L-5L-5L-5L-5L) and six 7-atom-layered slabs (7L-7L-7L-7L-7L-7L), leading to 30 atomic layers and 42 atomic layers respectively, have been investigated under strain. The calculated band structures are shown in Figure 13. As mentioned above, TIs exhibit gapless bands and Dirac-like dispersion, however a Dirac-like cone is only observed in PbBi 2 Te 4 ( Figure 13b). The Dirac-like cone observed in bulk Bi 2 Te 3 disappears in the Bi 2 Te 3 thin film (Figure 13a). The time-reversal invariant energy bands in 3D TI are characterized by 4 Z 2 invariants distinguishing 16 phases with 2 general classes: weak (WTI) and strong (STI) topological insulators [56,57]. Due to the spatial inversion symmetry of hexagonal Bi 2 Te 3 and PbBi 2 Te 4 , the determination of Z 2 invariants is simplified. Indeed, they can be obtained from the parity of the occupied states at the eight time-reversal invariant momenta Γ i in the Brillouin zone using the following formula [58]: where δ i = ∏ N m=1 ξ 2m (Γ i ) with N are the Kramers doublets in the valence band, and ξ 2m (Γ i ) is the parity of the 2m th occupied state at (Γ i ).

Conclusions
We have investigated the effect of various tensile and compressive strains on the electronic structure and transport properties of Bi 2 Te 3 , PbBi 2 Te 4 , PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 by means of first-principle calculations. Our calculations reveal that Bi 2 Te 3 , PbBi 2 Te 4 , PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 without strains are semiconductors. When applying strains, whether being tensile or compressive, a semiconductor-to-metal transition occurs for all the compounds, with the strain values corresponding to the transition being different for each compound. The range of non-zero band gaps is larger under compressive strains than under tensile ones. Seebeck coefficients' maxima were located at −2%, −6%, −3% and 0% strain for ptype Bi 2 Te 3 , PbBi 2 Te 4 , PbBi 4 Te 7 and Pb 2 Bi 2 Te 5 , respectively. These strain values correspond roughly to those related to the band-gap maxima for each compound. Besides, compressive strains induce larger Seebeck coefficients than tensile ones. Within the range of opengap, the electrical conductivity decreases as the compressive strain increases. Electrical conductivity is more sensitive to the applied strain for p-type doping than for n-type doping. Based on the analysis of charge density differences and the quantum theory of atoms in molecules, we have found that for the van der Waals interactions located between the slabs, the charge density increases and the bond degree goes from a positive value to a negative one when going from a compressive strain to a tensile one.
Furthermore, the topological properties of Bi 2 Te 3 (under −8.4% strain) and PbBi 2 Te 4 (under −2.2% strain) have been investigated. The strained Bi 2 Te 3 was found to be topologically trivial, whereas the strained PbBi 2 Te 4 can be estimated as a weak TI. Funding: The PhD thesis of W. Ma is financially supported by the China Scholarship Council (CSC). This work was granted access to the HPC resources of the Centre Informatique National de l'Enseignement Supérieur (CINES), Montpellier, France, under allocation A0070806881 made by the Grand Equipement National de Calcul Intensif (GENCI). It was also granted access to the HPC resources of Aix-Marseille University financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program "Investissements d'Avenir" supervised by the Agence Nationale de la Recherche.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: See supplementary data on MDPI website.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations were used in this manuscript: