Structural and physical properties of 99 complex bulk chalcogenides crystals using first-principles calculations

Chalcogenide semiconductors and glasses have many applications in the civil and military fields, especially in relation to their electronic, optical and mechanical properties for energy conversion and in enviormental materials. However, they are much less systemically studied and their fundamental physical properties for a large class chalcogenide semiconductors are rather scattered and incomplete. Here, we present a detailed study using well defined first-principles calculations on the electronic structure, interatomic bonding, optical, and mechanical properties for 99 bulk chalcogenides including thirteen of these crytals which have never been calculated. Due to their unique composition and structures, these 99 bulk chalcogenides are divided into two main groups. The first group contains 54 quaternary crystals with the structure composition (A2BCQ4) (A = Ag, Cu; B = Zn, Cd, Hg, Mg, Sr, Ba; C = Si, Ge, Sn; Q = S, Se, Te), while the second group contains scattered ternary and quaternary chalcogenide crystals with a more diverse composition (AxByCzQn) (A = Ag, Cu, Ba, Cs, Li, Tl, K, Lu, Sr; B = Zn, Cd, Hg, Al, Ga, In, P, As, La, Lu, Pb, Cu, Ag; C = Si, Ge, Sn, As, Sb, Bi, Zr, Hf, Ga, In; Q = S, Se, Te; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {x} = 1$$\end{document}x=1, 2, 3; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {y} = 0$$\end{document}y=0, 1, 2, 5; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {z} = 0$$\end{document}z=0, 1, 2 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {n} = 3$$\end{document}n=3, 4, 5, 6, 9). Moreover, the total bond order density (TBOD) is used as a single quantum mechanical metric to characterize the internal cohesion of these crystals enabling us to correlate them with the calculated properties, especially their mechanical properties. This work provides a very large database for bulk chalcogenides crucial for the future theoretical and experimental studies, opening opportunities for study the properties and potential application of a wide variety of chalcogenides.


Methods
In this work, two well-defined density functional theory (DFT) based methods were used: the Vienna Ab initio Simulation Package (VASP) 43,44 and the orthogonal linear combination of atomic orbitals (OLCAO) method 45 . VASP was used to optimize the crystal structures and to calculate the mechanical properties. The Perdew-Burke-Ernzerhof (PBE) generalized-gradient-approximation (GGA) 46 was used for the exchange and correlation potential in solving the Kohn-Sham equation with an energy cut-off of 500 eV. This energy cut-off has been carefully tested for all crystals and it is found to be the best choice balancing the fast convergence and accurate ground state energy. For the VASP calculation, the Monkhorst scheme 47 with different k-point mesh ranging from 9 × 9 × 5 and 5 × 9 × 5 was used for Ba 4 AgInS 6 (48 atoms) and Ag 2 In 2 SiSe 6 (44 atoms) respectively to 8 × 9 × 9 for medium sized crystals such as Cu 2 ZnSiSe 4 (16 atoms). The electronic and ionic force convergence criteria were set at 10 −6 eV and 10 −4 eV/Å respectively. The VASP optimized structures for all crystals were used as input to calculate the electronic structure, interatomic bonding, and optical properties using the OLCAO method with different choice of basis size for each atom in the database 45 . The use of localized atomic orbitals in the basis expansion in contrast to the plane-wave expansion is particularly effective for calculating the electronic structure and interatomic bonding for both crystalline [48][49][50][51][52] and non-crystalline materials 53,54 especially those with complex structures typical in the biomolecular systems 55,56 . In the OLCAO calculation, a sufficient number of k-points mesh were used for band structure and density of states (DOS) calculations based on the size of crystal. For interatomic bonding, a more localized minimal basis (MB) set is used based on Mulliken scheme 57 . Equations (1) and (2) show the formulae for effective charge(Q * α ) and bond order (BO) values, or the overlap population ρ αβ between any pair of atoms ( α, β).
In the above equation, S iα,jβ are the overlap integrals between the ith orbital in α th atom and the jth orbital in β th atom. C m jβ is the eigenvector coefficients of the mth occupied band. The BO from Eq. (2) defines the relative strength of the bond. The summation of all BO values in the crystal gives the total bond order (TBO). When normalized by the cell volume, we obtain the total bond order density (TBOD). TBOD is a single quantum mechanical metric to describe the internal cohesion of the crystal 58 . It can be conveniently decomposed into partial components or the partial bond order density (PBOD) for any structural units or groups of bonded atoms. The details for the calculation of interatomic bonding, optical, and mechanical properties are described in the Sects. 1 and 2 in the Supplementary Information (SI). The orthogonalized linear combination of atomic orbitals (OLCAO) method was used to calculate the electronic structures and partial charge distributions of the chalcogenide crystals under study. The OLCAO method is an all-electron method based on the local density approximation of DFT. It uses the atomic orbitals expanded in Gaussian-type orbitals (GTO) in the basis expansion. This method is particularly efficient for calculating the electronic structure of different systems. In the present calculation, a full basis set, which consisted of the core orbitals, occupied valence orbitals, and the next www.nature.com/scientificreports/ empty shell of unoccupied orbitals for each atom, was used for the self-consistent potential and the density of states (DOS) calculations. The more localized minimal basis was used for partial charges calculation under the Mulliken scheme. The partial charge Q on each atom is defined as the deviation of charge from the calculated effective charge Q * from the neutral atom charge ( Q 0 ) or Q = Q 0 − Q * . So negative δQ implies gain of electrons or electronegative and positive Q means loss of electron or electropositive. Table S1 in the SI lists the 99 quaternary crystals of the first group (A 2 BCQ 4 ) (1-54) and the second group (A x B y C z Q n ) . The data presented includes the name of the crystal, its space group, the VASP relaxed crystal parameters, and the available experimental parameters with appropriate references. In all subsequent discussions, the same specific order and the ID number labeled for the crystal are maintained to avoid any confusion. To make it easier for the readers to identify the specific crystal, we add the ID number in front of the crystal name most of the time. In both groups, we distinguish the chalcogen elements by referring them as Q, and the other elements by A, B, and C. Any exceptions will be specifically mentioned.

Results and discussion
Electronic structure. The most important quantity in a crystal to understand its physical and chemical properties is its electronic structure. The calculated results for the band structures and the density of states (DOS) for the 99 chalcogenide crystals are shown in Figs. S1 and S2 respectively. Partial charge. An important electronic structure property in a crystal is the distribution of the partial charge (PC) on each atom in the crystal. PC is the deviation of the effective charge Q * from the neutral atom charge Q 0 or simply the charge transfer. The calculated Q * on each atom in the 99 chalcogenide crystals is listed in column 3 of Table S3. The generally accepted concept is that chalcogen elements (S, Se, Te) receive electrons from the other atoms in the chalcogenide compounds or they are electronegative. The other atoms in a chalcogenide compound (except for Cu and Ag in most cases) lose charge and they are mostly electropositive. It turns out that this generally accepted notion is grossly over-simplified, and the real situation is far more complicated and subtle. With 99 chalcogenide crystals formed by 30 elements, we are in the unique position to analyze the PC distribution and charge transfer mechanism in much greater detail than for just one or a few crystals. In the first group with the structure (A 2 BCQ 4 ), the PC values for most of A elements are negative similar to the chalcogen elements except for some crystals such as 36-Ag 2 ZnSiS 4 , 38-Ag 2 ZnSnSe 4 , and 39-Ag 2 CdGeS 4 where the A element (Ag) loses charge. While the B and C elements have positive PC. To explain why Cu and Ag are electronegative in most of the crystals in group one, we must dig deeper. Cu and Ag are noble metals with full shell of 10 3d and 10 4d electrons respectively, and metals usually lose charge to the other elements in the chemical compound, specially to the nonmetals elements, but our calculations for the partial charge and effective charge indicate that these two elements (Cu, Ag) can gain or loss charge depending on chalcogen elements their interact. First, for crystals that contain Cu and Ag with the nonmetal chalcogen elements (S, Se) and later with Te (a metalloid element) which is very different than S and Se. The bonding properties indicate that Cu and Ag have two types of bonds. depending on the chalcogen elements (S or Se). They form bonds that are ionic bonds (metal-nonmetal bond), and between two metal elements (Cu, Ag). Cu and Ag lose charge to the nonmetallic chalcogen elements S and Se and these bonds occur at shorter bond lengths (1.5-2.0 Å). The other type of bonds is between Cu or Ag and with one of the metal elements Zn and Cd and with one of the alkaline earth elements (Mg, Sr, Ba). These bonds occur at long bond lengths ( > 3.8 Å www.nature.com/scientificreports/ bar in each column indicate the averaged PC for that element. Figure S3 shows the PC distribution of the first 13 elements. Depending on the crystals they originate, they generally have a wide range of scatted values and only Cu has an average negative PC of − 0.104 e − and Ag has a slightly positive average charge of 0.026 e − . From Zn to Pb, the average PC range from 0.310 e − in Zn to 0.811 e − in Cd with most of them having close average values. An important facet is a wide scattering of PC values for individual crystals indicating the averaged PC value for any element in chalcogenide crystals has no real meaning since it depends on the crystal components and the overall charge transfers between them. The situation is slightly different in Fig. S4 for the 14 elements from P to Lu mainly because there are much less crystals involving these group of elements. One of the spectacular outliers of the rare earth element La from a single crystal 74-Ba 2 LaGaSe 5 . Other outliers include Ba which consists of 11 crystal distinctively separated into 2 groups with average PC of 1.729 e − between them. The most important PC information is the distribution of the three chalcogen elements in Interatomic bonding. An important characteristics is to investigate the interatomic bonding between every pair of atoms in the crystal represented by the BO value that quantifies the strength of the bond 58 . These data can then be used to obtain the key parameter TBOD discussed in the method section. The highest TBOD among these 99 crystals is in 1-Cu 2 ZnSiS 4 followed by 32-Cu 2 MgSiS 4 , while the lowest TBOD is from 77-Ba 4 A-gInS 6 closely followed by 72-Ba 2 GaBiSe 5 (see Table S3). In Figs. 5 and 6, we display the TBOD for the 99 crystals in 2 main groups respectively. The first group with 54 crystals (A 2 BCQ 4 ), is divided into two subgroups (Cu- www.nature.com/scientificreports/ related and Ag-related) indicated with horizontal bars in Fig. 5, while the second group (A x B y C z Q n ) with 45 crystals in Fig. 6 is divided into 9 subgroups indicated by horizontal bars which are related to the elements A= (Cu, Ag, Ba, Cs, K, Li, Lu, Sr, Tl) respectively. It can be seen that the first main group A 2 BCQ 4 have relatively higher TBOD values than the second main group (A x B y C z Q n ). In Fig. 5 for the first main group, crystals with A = Cu, Ag; B = Hg; C = Sn; Q = Se, Te have the lowest TBOD. The data in Fig. 6 for the second group are more scattered. However, we can still observe that the Cu-related crystals have the highest TBOD whereas the Ba-related crystals have the lowest TBOD. The Ag-related, Cs-related, and Li-related crystals all have very close TBOD values. An important observation from the TBOD calculations is that in moving from S to Se to Te through these 99 crystals, the TBOD value decreases. The possible reason has to do with the size of the atoms which affects the volume of the crystal resulting in longer bond lengths hence smaller TBOD. That being said, we cannot say that the chalcogen elements completely control the TBOD value of the crystal because the elements A, B and to a lesser extend C also play a role in determining their total bond order and the volume of the crystal, hence the www.nature.com/scientificreports/ TBOD. We need to reiterate that the use of the novel concept of TBOD is extremely helpful to extract the internal cohesion of otherwise very complex crystals of chalcogenide compounds.
Optical properties. The optical properties of these 99 chalcogenide crystals are particularly important because of their many optical applications. They can be calculated relatively easily within the one-electron random phase approximation using the OLCAO method. It starts with the calculations of the imaginary part of the www.nature.com/scientificreports/ dielectric function ǫ 2 (ω) from the interband optical transition with explicit inclusion of the dipole matrix elements. The real part of the dielectric function ǫ 1 (ω) is obtained from ǫ 2 (ω) by applying Kramers-Kronig transformation. The refractive index n of the crystal is obtained as the square root of ǫ 1 (ω) in the zero-frequency limit, or n = √ ǫ 1 (0) . More details can be found in Sect. 1 in SI including the calculation of the energy loss function (ELF) and the evaluation of the Plasmon frequency ω p . The calculated real and imaginary parts of the dielectric functions for these 99 crystals are shown in Fig. S6, and the refractive indices n are listed in Table S4. As can be seen, each crystal has its unique absorption features or interesting optical properties. We selectively discuss some of them in A 2 BCQ 4 (with A = Ag, Cu; B = Sr, Ba; C = Ge, Sn; Q = S, Se) and also the following Ba-related and Sr related crystals (71-Ba 2 AlSbS 5 , 72-Ba 2 GaBiSe 5 , 73-Ba 2 AsGaSe 5 , 74-Ba 2 LaGaSe 5 , 75-Ba 2 LuGaSe 5 , 76-Ba 2 LuInSe 5 , 77-Ba 4 AgInS 6 , 86-SrCdGeS 4 , 87-SrCdGeSe 4 ). They are part of the results shown in Fig. S6. The general feature for these selected crystals is that ǫ 1 (ω) reaches a peak at certain energy and then starts decreasing, eventually reaching a negative value. It then increases again before finally leveling off to zero value again. Since the crossing at zero value is generally associated with Plasma excitation, this implies that in these group of selected crystal, they may have two plasma excitations and two Plasma frequencies, such as that reported recently in Ref. 59 . This unique feature is not specifically related to the structure of the crystals (A 2 BCQ 4 or A x B y C z Q n ), but strongly related to the presence of Ba and Sr in these crystals. This is an important observation revealed in the present study which clearly originate from the electronic structure. From the dielectric function, we can obtain the ELF for all the 99 chalcogenides which are shown in Fig. S7. The ELF is an important part of the optical properties of crystals and represents the collective excitation of excited electrons at high frequency. The position of its main peak is defined as the Plasmon frequency ω p which usually occurs at the frequency when the real part of dielectric function vanishes. Below ω p , the incident waves will be mostly reflected. The ω p values for the 99   www.nature.com/scientificreports/ crystals are listed in Table S3. The highest ω p is identified to be in 85-Lu 5 GaS 9 at 31.10 eV, while the lowest ω p is in 63-Ag 2 SnSe 3 at 15.50 eV.

Mechanical properties.
Mechanical properties of chalcogenide materials are much less studied compared to optical properties even though they are critical in many practical applications 60,61 . In this work, we fill this void by calculating the mechanical properties of the 99 chalcogenides using the method described in Sect. 2 in SI. The calculated elastic constants are listed in Table S5 and the mechanical parameters derived from them are listed in Table S6. The principal coefficients C 11 , C 22 , and C 33 , related to the unidirectional compression along the principal x, y and z directions, reflect the isotropic elasticity of the crystals and they are usually close in cubic crystals with high symmetry 62 . Equivalently, we say C 11 , C 22 , and C 33 reflect the resistance to the deformation of the crystal along x, y, and z directions. The elastic constants depend on the crystals structure (lattice parameters) and the strength of the bonds between the elements in the crystal. Table S5 shows that some of the crystals have C 22 higher than C 11 and C 33 which simply suggest that they are more compressible along x-and z-axes than along y-axis. We note that the two crystals 38-Ag 2 ZnSiS 4 and 63-Ag 2 SnSe 3 have negative elastic constants which implies a negative stiffness. This behavior could be related to a material instability 63 . To the best of our knowledge, this is the first time that this observation is reported for these two crystals either experimentally or computationally. From the elastic coefficient Cij and the compliance tensor Sij (Sij = 1/Cij), the mechanical parameters bulk modulus (K), shear modulus (G), Young's modulus (E), and Poisson's ratio ( η ) can be obtained using Voigt-Reuss-Hill (VRH) approximation for polycrystals as explained in SI [64][65][66] . They are listed in Table S6. Another useful parameter is G/K or Pugh's modulus ratio 67,68 . According to Pugh's criterion, crystals with G/K larger than 0.57 are brittle, and those with less than 0.57 tend to be more ductile [68][69][70][71] . Figure 7 shows the scattered plot of G versus K for the 54 crystals in the first main group, while     We also notice that K and G values decrease in moving from S to Se to Te. This feature can again be attributed to the larger size of Te atom resulting in weaker bond strength and more brittle crystals. The same feature has been observed for the TBOD and this implies a strong correlation between the mechanical properties and TBOD. Figures 13 and 14 show the scattered plots of TBOD versus Poisson's ratio η to explore any trend for these 99 crystals in two groups. It appears that the relationship between TBOD and Poisson's ratio is not that much easy to explain and the data are scattered. In a broader scale, we can still observe a general trend that a larger TBOD implies a larger Poisson's ratio. Identification of the underlying correlation between the electronic structure and the bonding characteristics of the chalcogenide crystals is one of the main objectives in this work. In this regard, exploring the connection between mechanical properties of these 99 crystals and the TBOD could be revealing. In Figs. 15,16,17, and 18, we plot the bulk modulus and shear modulus versus TBOD for the 99 crystals in the two separate groups. For the first main group of 54 crystals, there is a very clear correlation between the two moduli (K, G) and the TBOD. K and G increase linearly with the TBOD with only a few outliers. For the second main group of 45 crystals, the situation is slightly different. Moreover, we divided the crystals in this main group into two subgroups. In the first subgroup (71-Ba 2 AlSbS 5 , 72-Ba 2 GaBiSe 5 , 73-Ba 2 AsGaSe 5 , 74-Ba 2 LaGaSe 5 , 75-Ba 2 LuGaSe 5 , 77-Ba 4 AgInS 6 , 91-Tl 2 PbZrS 4 , 92-Tl 2 PbZrSe 4 , 93-Tl 2 PbHfS 4 , 94-Tl 2 PbHfSe 4 -mostly contains the Ba-related and some of the Tlrelated crystals inside the blue circle), they are quite scattered with no clear correlation for K and G with TBOD. While in the second subgroup (the remaining crystals), there is a reasonable correlation with K and G steadily increasing linearly with TBOD. The crystal 85-Lu 5 GaS 9 deviates from this behavior. It exhibits a sharp increase of G and K with TBOD. We conclude that the mechanical properties of these 99 chalcogenides are strongly related to the strength of the interatomic bonds which are collectively represented by the single matric TBOD.

Conclusion
In conclusion, a comprehensive library of the electronic structure, interatomic bonding, optical, and mechanical properties of the 99 chalcogenides was assembled on the basis of extensive first-principles calculations. We summarize below the new insights obtained and the conclusions reached in this study by using the same approach and resulting in a pronounced consistency of the result.       www.nature.com/scientificreports/ TBOD. In the second main group, the Cu, Cs-related crystals have the highest TBOD and the Ba-related crystals have the lowest TBOD. The Ag-related, Cs-related and Li-related crystals have very similar TBOD. 5. The optical spectra of the 99 crystals are calculated within the random phase approximation. They have medium values for the refractive index n in the range of 4.424-4.524. 6. The mechanical properties for these 99 crystals show large variations between the two main groups (A 2 BCQ 4 ) and (A x B y C z Q n ). The first group crystals tend to be more ductile than those from the second group. They are both reasonably correlated with the TBOD. The results and correlations summarized above can facilitate the design of new chalcogenide crystals and glasses for potential novel applications. They constitute a formidable database for complex chalcogenide crystals that was not available before. We believe that the present work can facilitate the discovery and production of good-quality chalcogenides with a variety of energy or environmental applications.

Data availability
All the data in this paper including those in the supplementary materials are freely available by contacting the one of the corresponding authors (chingw@umkc.edu) Received: 21 November 2020; Accepted: 6 April 2021