Generative design of stable semiconductor materials using deep learning and density functional theory

Abstract

Semiconductor device technology has greatly developed in complexity since discovering the bipolar transistor. In this work, we developed a computational pipeline to discover stable semiconductors by combining generative adversarial networks (GAN), classifiers, and high-throughput first-principles calculations. We used CubicGAN, a GAN-based algorithm for generating cubic materials and developed a classifier to screen the semiconductors and studied their stability using first principles. We found 12 stable AA${}^{\prime}$MH₆ semiconductors in the F-43m space group including BaNaRhH₆, BaSrZnH₆, BaCsAlH₆, SrTlIrH₆, KNaNiH₆, NaYRuH₆, CsKSiH₆, CaScMnH₆, YZnMnH₆, NaZrMnH₆, AgZrMnH₆, and ScZnMnH₆. Previous research reported that five AA${}^{\prime}$IrH6 semiconductors with the same space group were synthesized. Our research shows that AA${}^{\prime}$MnH₆ and NaYRuH₆ semiconductors have considerably different properties compared to the rest of the AA${}^{\prime}$MH₆ semiconductors. Based on the accurate hybrid functional calculations, AA${}^{\prime}$MH₆ semiconductors are found to be wide-bandgap semiconductors. Moreover, BaSrZnH₆ and KNaNiH₆ are direct-bandgap semiconductors, whereas others exhibit indirect bandgaps.

Introduction

Semiconductors are essential components of modern devices that use transistors, light-emitting diodes¹, integrated circuits², photovoltaic³, solar cells⁴, and so on^5,6,7. Semiconductors exhibit variable resistance since electron flow can be controlled by light and heat. Therefore, these materials can be used for energy conversion, and digital switching⁸. The elemental semiconductors found from Group XIV in the periodic table, like Si and Ge, and the compounds of Ge are widely used in electronics, photovoltaic and optoelectronic devices. However, semiconductors with various properties are required for industrial applications^8,9. For instance, good thermal conductivity and electric field breakdown strength, and also wide bandgap of SiC semiconductor make it a suitable material for high-temperature, high-power, high-frequency, and high-radiation conditions¹⁰. Thus, computational approaches for exploring semiconductors are essential to enhance future technologies. High-throughput screening with the aid of first-principles calculations was performed by several groups to discover optoelectronic semiconductors. Setyawan et al. and Ortiz et al. reported the high-throughput screening and data-mining frameworks to investigate bandgap materials for radiation detection^11,12,13. High throughput material screening by Zhao et al. found that Cu-In-based Halide Perovskite as potential photovoltaic solar absorbers^13,14. Based on 4507 hypothetical materials, Li et al. suggest 23 candidates for light-emitting applications, and 13 potential compounds for solar cell technologies^13,15. Such examples indicate that high-throughput screening can now be used to explore promising semiconductor materials.

Generative adversarial networks (GANs) are a kind of generative models that learn patterns/distribution from input data¹⁶. GANs use two sub-models to train a generative model. The generator model generates fake data, and the discriminator model learns to tell fake data from real data. The two sub-models are trained simultaneously to achieve a Nash Equilibrium: the generator can generate data that the discriminator can recognize half the chance. Wasserstein distance¹⁷ and gradient penalty¹⁸ are introduced during training in order to overcome mode collapse and improve the training stability in original GANs¹⁶. There are a limited number of works that leverage GANs to generate crystal structures in material science. The reasons behind that are: 1) Crystal structures have so many formations, such as a different number of elements and number of atoms in a unit cell. It is hard to come up with a unified representation to make GANs learn from them like images or text; 2) GANs used in computer vision cannot generate crystal structures that satisfy physics or symmetric constraints. For instance, GANs easily generate materials that are not recognizable or that have crowd atoms in a unit cell. CrystalGAN¹⁹ is believed to be the first work that uses GANs to generate materials. It applies CyClyGAN²⁰ to simple systems mapping ternary a hydride into another. In²¹, Kim et al. use WGAN-GP¹⁸ to train a generative model to generate Mg-Mn-O systems with atom coordinates as the input. All the works above only consider a simple or specific family of materials at a limited scale. CubicGAN proposed by Zhao et al.²², however, is the first work that generates materials at a large scale.

In this research, we developed a binary classifier to filter the semiconductors/Insulators (nonmetals) from the dynamically stable quaternary Cubic materials discovered using the CubicGAN model, where high-throughput calculations were done with the assistance of a GAN model and density functional theory (DFT). We studied the most important elemental and electronic properties, which are helpful to distinguish the nonmetals and metals using the machine learning models. In addition, we carried out DFT calculations for those semiconductors to corroborate the thermodynamic stability and semiconductor properties. As a result, we find that 12 cubic semiconductors of a particular class of materials, which we label as AA${}^{\prime}$MH₆, are thermodynamically stable against their competing phases. We further performed the DFT calculations to study their structural, mechanical, thermodynamic, and electronic properties. Our results show that AA${}^{\prime}$MnH₆ and NaYRuH₆ have higher C_ii (i = 1, 2, 3) elastic constants, bulk modulus, shear modulus, and Young’s modulus compared to the respective mechanical properties of the rest of the AA${}^{\prime}$MH₆ materials. At temperatures less than 200 K, AA${}^{\prime}$MnH₆ and NaYRuH₆ have lower specific thermal capacity (C_v) relative to other AA${}^{\prime}$MH₆ materials. The highest C_v at 300 K found in this work is from BaSrZnH₆ (127.96 JK⁻¹mol⁻¹). Moreover, hybrid functional calculations show that all AA${}^{\prime}$MH₆ materials are wide-bandgap semiconductors, which will be useful to develop optical and high-temperature power devices^23,24.

Results and discussion

Dataset of nonmetals and metals

As the CubicGAN model generates only ternary and quaternary materials, we first analyzed the number of nonmetals (semiconductors and insulators), and metals in the material project (MP) database²⁵, as shown in Table 1. We collected all the ternary and quaternary materials, where the bandgap details are available, using the Pymatgen code²⁶. It could be found that ≈44 % of the ternary materials are nonzero bandgap materials while ≈ 56 % are metals. However, ≈73 % of the quaternary materials are semiconductors or insulators, whereas only ≈27 % of them are metals. This indicates that the probability of finding a stable quaternary material with a nonzero bandgap is higher compared to finding that in a ternary material set. We also compared the same details of the cubic materials. Interestingly, ≈80 % of the cubic ternary materials are metals, and only ≈20 % of them are nonmetals. On the contrary, the quaternary cubic materials have 30 % more nonzero bandgap materials than the number of metals. It shows that there is a low probability of discovering a nonzero bandgap cubic ternary compound. Instead, in this project, we mainly focused on the quaternary cubic materials for finding stable semiconductors. In this way, by reducing the search space of the materials, we can shorten the computational time taken by the DFT calculations.

Table 1 The total number (N) and percentage of ternary and quaternary nonmetals (NM) and metals (M) for all the types of Bravais lattice and also for Cubic structures (Cubic-NM and Cubic-M) in the Material project database²⁵.

Feature importance

Understanding which features are significant during the classification will be vital for discovering semiconductors. In Section 2.1, we could show that quaternary materials have a higher percentage of semiconductors compared to the ternary materials. Next, we analyzed which features have higher importance than others for classifying a quaternary material as metal or nonmetal. Feature importance (FI) of random forest algorithm is defined as the mean of the impurity decrease within each tree. This built-in feature of the random forest makes it convenient and a widely used method to calculate FI. Here, we trained our RFC model for the whole quaternary materials data set. The classification report of this model is in Supplementary Information. Even though both Avg. and the maximum difference of each atomic/electronic property were considered for the RFC model, only three features related to maximum difference have FI greater than 1 %. This indicates that Avg. value of the properties plays a significant role when classifying a material as metal or nonmetal. The top features of FI ⪆ 2.0% are mentioned in Fig. 1. Avg. Availability of metallic elements has the highest FI, while Avg. availability of nonmetal also has a FI of around 2 %. This indicates that having a metallic or nonmetallic element is important for the material to be a metal or a semiconductor/insulator. It is generally accepted that metallic elements have a higher boiling point and higher density compared to that of nonmetals. It should be noted that the elemental properties like metallicity, being semiconductor/insulator, density, and boiling point are properties of the bulk material formed with a given element. Since the availability of metallic and nonmetallic elements plays a significant role, the boiling point and density of those elements also can become important features when classifying metals and nonmetals. It is also clear that electronic properties like Avg. number of unfilled orbitals, Avg. number of p-valence electrons, and Avg. availability of +2 and +3 oxidation states have high FI.

Fig. 1: Feature importance as a percentage ( ⪆ 2.0%) from the RFC model.

Labels on the x-axis: metal atom (M), +2 oxidation state (+2), boiling point (T_boil), atomic density (ρ), number of p-valence electrons (Np_val), +3 oxidation state (+3), nonmetal atom (NM), number of unfilled valence orbitals (N_unfilled). Avg. and Avail. stands for weighted average and availability, respectively.

We also studied the descriptors to understand how the number of metals and nonmetals depends on the percentage availability of the metal (P_M), nonmetal (P_NM) and transition-metal (P_TM) elements in the chemical formulas. We use M, NM, and TM to indicate the type of elements to avoid confusion between material class (metal or nonmetal) and element type (metal, nonmetal, transition metal). Figure 2 shows the violin plots with all the 39024 quaternary materials against those three atomic properties. Here, P_M = 100%, P_NM = 100%, and P_TM = 100% for a given chemical formula when all the elements are M, NM, and TM, respectively. Figure 2(a) clearly evidences that nonmetals dominate until P_M ≈ 60 %. The ratio between amounts of metals and nonmetals (metals : nonmetals) is around 1: 3 at P_M < 60 %. This becomes approximately 5: 1 after 60 %, showing the probability of finding a semiconductor/insulator decreases. On the contrary, Fig. 2(b) shows the opposite behavior of metals and nonmetals, while P_NM alters. Moreover, it is clear that semiconductors and insulators prefer a lower number of TM elements relative to the other element types. At P_TM > 30 %, number of metals become significant compared to that of nonmetals. When P_TM ≤ 5 %, metals : nonmetals ratio is 1: 6.

Fig. 2: Number of materials in each material class (metal or nonmetal) as a function of percentage availability (P) of element type in the chemical formula.

Here, we considered (a) metal:M, b nonmetal:NM, and c transition-metal:TM elements.

Predicting Semiconductors

We further analyzed the error of the DNN and RFC models trained with quaternary cubic materials data. The 10-fold cross-validation accuracy results for each training step of the DNN model are 0.86, 0.92, 0.91, 0.97, 0.94, 0.94, 0.88, 0.94, 0.94, 0.86. Those of the RFC model are 0.86, 0.89, 0.88, 0.9, 0.87, 0.88, 0.90, 0.90, and 0.88. Thus, the mean accuracy was obtained for the DNN (RFC) model as 0.92 ± 0.034 (0.88 ± 0.013). Figure 3 shows the normalized confusion matrices for the classifiers. It is apparent that 33 (32) % of the instances were classified as true metals while 65 (60) % of the materials were listed as true nonmetals by the DNN (RFC) classifier. The percentages of false metals and false nonmetals from the DNN (RFC) model were 9.8 (4.9) % and 1.2 (2.5) %, respectively. The classification report for the model is shown in Table 2. It is clear that the DNN (RFC) classifier predicts whether a quaternary material is a metal or nonmetal with 0.88 (0.91) accuracy. Precision is the matrix that compares the number of true positive instances with the number of predicted positive instances. In our work, the DNN (RFC) model classifies a material as a nonmetal with 0.76 (0.96) and metal with 0.76 (0.84) precision. The recall is a measure of the number of correctly predicted positive cases compared to the total number of positive cases in the dataset. Table 2 shows that there is 0.85 (0.91) recall for nonmetal, while there is 0.93 (0.93) recall for metals from the DNN (RFC) model. By combining precision and recall, F1-score can be calculated as 0.90 (0.93) for nonmetal and 0.84 (0.88) for metal classes. Therefore, the predictions of semiconductors/insulators from our DNN and RFC models can be expected to be highly accurate.

Fig. 3: Normalized confusion matrices of the classifiers for metals (M) and nonmetals (NM).

The confusion matrix for deep neural network (DNN) is in purple, while that of random forest (RFC) is in blue.

Table 2 The classification report of the test set for the nonmetal/metal classifiers.

As seen in Table 2, the RFC model exhibits a slight improvement over the DNN model. To show the methodology of finding stable semiconductors based on generative adversarial networks, we applied our RFC classifier on CubicGAN predicted mechanically and dynamically stable quaternary materials. Out of 323 quaternary materials predicted by the CubicGAN model, 137 compounds were classified as nonmetals.

Structure and thermodynamic stability

We carried out our DFT calculations on those nonmetals to find thermodynamically stable semiconductors. We discovered that 12 semiconductors, which have chemical formulas in the form of AA${}^{\prime}$MH₆, exhibit zero energy-above-hull against the respective competing phases. Those are BaNaRhH₆, BaSrZnH₆, BaCsAlH₆, SrTlIrH₆, KNaNiH₆, NaYRuH₆, CsKSiH₆, CaScMnH₆, YZnMnH₆, NaZrMnH₆, AgZrMnH₆, and ScZnMnH₆. We also find that Kadir et al. reported 5 different AA${}^{\prime}$MH₆ type semiconductors, where M = Ir²⁷. They were able to synthesize NaCaIrH₆, NaSrIrH₆, NaBaIrH₆, KSrIrH₆, and KBaIrH₆ by direct combination of the alkali (Na and K), alkaline earth (Ca, Ba, and Sr) binary hydrides/deuterides with Ir powder. Their X-ray and neutron powder diffraction studies confirm that those semiconductors have the space group symmetry F-43m. Furthermore, the open quantum materials database (OQMD)^28,29 contains the structural properties and band gaps of NaCaIrH₆, NaSrIrH₆, NaBaIrH₆ semiconductors and the MP database has those information on NaCaIrH₆, and NaBaIrH₆ semiconductors²⁵ (See Supplementary Information).

CubicGAN generates conventional structures with cubic Bravais lattice with F-43m (216) space group for AA${}^{\prime}$MH₆ materials, which have 36 atoms. On the contrary, the primitive unit cell with hexagonal Bravais lattice has only 9 atoms. Therefore, we considered the hexagonal unit cell to lower the computational time of the DFT calculations. In the primitive unit cells (see Fig. 4), green and red sites are symmetrically equivalent, while grey sites are located in the right middle of the hexagonal unit cell. Thus, we label the green and red sites as A and A${}^{\prime}$, while the middle site is M. Rest of the 6 sites are occupied by H atoms. In the research work of Kadir et al., they considered alkali atoms as A atoms, alkaline earth atoms as A${}^{\prime}$ atoms, and M atoms as Ir. In this research, our findings show that both A and A${}^{\prime}$ atoms can be alkali atoms (E.g., CsKSiH₆) or alkaline earth atoms (E.g., BaSrZnH₆). Moreover, the M atom can be a transition metal atom or even Al or Si. Therefore, our experiments show that those materials can have high chemical diversity.

Fig. 4: Side view of the structure of AA${}^{\prime}$MH₆ materials with primitive and conventional unitcells, which are indicated by red lines.

The green, purple, grey and blue spheres indicate the A, A’, M and H atoms, respectively.

Fig. 5: Phonon dispersion of AA${}^{\prime}$MH₆ materials.

The positive phonon frequencies in the figure indicate that the materials are dynamically stable.

The lattice parameters, A-H, M-H, A-M, and A-A${}^{\prime}$ bond lengths, are mentioned in Table 3. The primitive hexagonal unit cells have a/c = 1 lattice parameter ratio making a = b = c. As shown in Table 3, Mn-related AA${}^{\prime}$MH₆ and NaYRuH₆ structures have the shortest lattice parameters compared to the rest of the materials. They have lattice parameters less than 5.0 Å, while other materials have greater than 5.4 Å. All A, A${}^{\prime}$ and M elements make bonds with H atoms. A and A${}^{\prime}$ elements are bonded to twelve equivalent H atoms to form AH₁₂ and A${}^{\prime}$H₁₂ cuboctahedra. And also, M atoms make MH₆ octahedra by making bonds with 6 H atoms. An AH₁₂ (A${}^{\prime}$H₁₂) cuboctahedra shares corners with twelve equivalent AH₁₂ (A${}^{\prime}$H₁₂) cuboctahedra. Moreover, they share faces with four MH₆ octahedra³⁰. Due to symmetry, A-H and A${}^{\prime}$-H bond lengths are equal. M-H bond lengths are the shortest compared to other bonds for a given compound. A-A${}^{\prime}$ of Mn-related AA${}^{\prime}$MH₆ and NaYRuH₆ structures are less than 3.4 Å, and A-M and A${}^{\prime}$-M distances are less than 3.1 Å. It can cause strong interactions between those atoms. A-A${}^{\prime}$ distance for the rest of the materials is greater than 3.8 Å, and A-M and A${}^{\prime}$-M distances are greater than 3.3 Å, indicating relatively weaker interactions.

Table 3 The lattice parameter (a) and bond lengths in Å, and the formation energies in eV/atom for AA${}^{\prime}$MH₆ materials.

The thermodynamic stability of the AA${}^{\prime}$MH₆ materials against their elements was studied using the formation energies, which were based on the following equation.

$${E}_{{{{\rm{form}}}}}=\frac{1}{N}\left({E}_{{{{\rm{tot}}}}}-\mathop{\sum}\limits_{i}{x}_{i}{E}_{i}\right)$$

(1)

Here, E_tot is the total energy per unit formula of the material. x_i is the number of atoms of each element in the unit formula; i.e., 1 for A, A${}^{\prime}$, M atoms and 6 for H. N = ∑x_i; i.e., 9 for AA${}^{\prime}$MH₆. To find the atomic energies (E_i), we collected the most stable structures of each element using the Pymatgen code²⁶. Same DFT settings were used to calculate the energy of each element. It is clear that all the six materials have negative formation energies, which confirms their stability. We also carried out spin-polarized calculations for the AA${}^{\prime}$MH₆ semiconductors with transition metal atoms to reveal whether they form magnetism. We observed that those materials do not have magnetic groundstates. Thus, all the AA${}^{\prime}$MH₆ semiconductors are nonmagnetic materials.

Mechanical properties and stability

Next, we studied the mechanical properties and stability of the AA${}^{\prime}$MH₆ materials by calculating the elastic constants using the DFPT method. To analyse the mechanical properties, we used the Vaspkit code³¹, which computes the elastic constants by considering the AA${}^{\prime}$MH₆ cubic system. Since cubic unitcells has a = b = c lattice lengths and α = β = γ = 90⁰ lattice angles, C₁₁ = C₂₂ = C₃₃, C₄₄ = C₅₅ = C₆₆, and C₁₂ = C₁₃ = C₂₃³². Therefore, we mention only the three independent elastic constants (C₁₁, C₁₂ and C₄₄) in Table 4. It is clear that AA${}^{\prime}$MH₆ materials have relatively higher C₁₁ for AA${}^{\prime}$MnH₆ and NaYRuH₆, compared to the other four materials in Table 4. As discussed before, the lattice constants and A-A${}^{\prime}$ bond lengths of AA${}^{\prime}$MnH₆ and NaYRuH₆ structures are considerably lower than that of the rest of the materials. As illustrated by Fig. 4, A-A${}^{\prime}$ bonds are aligned in a, b and c directions. C₁₁, C₂₂, and C₃₃ are parallel to the a, b and c directions, respectively. Therefore, higher C_ii (i = 1, 2 and 3) can be mainly due to the strong interactions between the A and A${}^{\prime}$ atoms. Born stability criteria for the cubic systems are C₁₁ − C₁₂ > 0, C₁₁ + 2C₁₂ > 0 and C₄₄ > 0³². It is clear from Table 4 that all the eight materials comply with the above requirements.

Table 4 The mechanical properties of the AA${}^{\prime}$MH₆ materials.

We also calculated the Bulk modulus (K), Young’s modulus (Y), and isotropic Poisson’s ratio (μ) based on the Hill approximation³³ as mentioned in Table 4. The smallest K values were found from CsKSiH₆ (16.615 GPa), while the largest value was calculated from AgZrMnH₆ (120.755 GPa). SrTlIrH₆ (21.915 GPa) provides the lowest Y, while NaZrMnH₆ (156.876 GPa) exhibits the maximum Y. It is clear that NaYRuH₆ and all the Mn-based materials have significantly larger K and Y values than that of the other six materials. This can be mainly because of high C_ii (i = 1, 2, and 3) formed due to strong A-A${}^{\prime}$ bonds. Because of low Y, NaYRuH₆ and Mn-based AA${}^{\prime}$MH₆ materials can be considered stiffer materials relative to the other six semiconductors. And also, they exhibit more resistance to compression due to high K. All the μ values of the AA${}^{\prime}$MH₆ materials are between 0.2 and 0.4. maximum μ was found from SrTlIrH₆. Thus, SrTlIrH₆ has considerably low Y and high μ. This indicates that SrTlIrH₆ semiconductor is less stiff due to small Y and more deformable elastically at small strains due to large μ.

Thermodynamic properties and dynamical stability

The temperature of the highest normal mode of a crystal is known as the Debye temperature θ_D. This can be obtained by employing Debye sound velocity (ν_D) as explained by Eq. (2). Debye sound velocity can be calculated using the longitudinal and transverse sound velocities, which can be determined based on K and G as shown in Eq. (4)³⁴. Here, N, V₀, and ρ are the number of atoms, volume, and density of the unicell, respectively. And also, h is Plank’s constant, and k_B is Boltzmann’s constant.

$${\theta }_{{{{\rm{D}}}}}=\frac{h}{{k}_{{{{\rm{B}}}}}}{\left(\frac{3N}{4\pi {V}_{0}}\right)}^{\frac{1}{3}}{\nu }_{{{{\rm{D}}}}}$$

(2)

$${\nu }_{{{{\rm{D}}}}}={\left[\frac{1}{3}\left(\frac{2}{{\nu }_{{{{\rm{l}}}}}^{3}}+\frac{1}{{\nu }_{{{{\rm{t}}}}}^{3}}\right)\right]}^{-\frac{1}{3}}$$

(3)

$${\nu }_{{{{\rm{l}}}}}={\left(\frac{3K+4G}{3\rho }\right)}^{\frac{1}{2}}\,\,{{{\rm{and}}}}\,\,{\nu }_{{{{\rm{t}}}}}={\left(\frac{G}{\rho }\right)}^{\frac{1}{2}}$$

(4)

Table 5 shows the respective ρ, ν_l, ν_t, ν_D and θ_D values for AA${}^{\prime}$MH₆ crystals. Debye temperature of NaYRuH₆ and Mn-based AA${}^{\prime}$MH₆ materials are significantly higher than that of other AA${}^{\prime}$MH₆ materials. As θ_D depends on K and G (see Eq. (4) and (2)), enhanced θ_D is due to the high K and G of those semiconductors.

Table 5 The density (ρ), longitudinal (ν_l), transverse (ν_t) and average (ν_D) sound velocity, Debye temperature (θ_D) and specific thermal capacity at 300 K (${C}_{{{{\rm{v}}}}}^{{{300\,{\rm{K}}}}}$) for the AA${}^{\prime}$MH₆ materials.

We also plotted C_v as a function of temperature T using the Phonopy code³⁵. C_v can be determined based on the following expression,

$${C}_{{{{\rm{v}}}}}=\mathop{\sum}\limits_{{{{\bf{q}}}}j}{k}_{{{{\rm{B}}}}}{\left(\frac{\hslash {\omega }_{{{{\bf{q}}}}j}}{{k}_{{{{\rm{B}}}}}T}\right)}^{2}\frac{\exp (\hslash {\omega }_{{{{\bf{q}}}}j}/{k}_{{{{\rm{B}}}}}T)}{{[\exp (\hslash {\omega }_{{{{\bf{q}}}}j}/{k}_{{{{\rm{B}}}}}T)-1]}^{2}},$$

(5)

where ω_qj is the phonon frequency for q wave vector at jth phonon band index and ℏ is the reduced Plank’s constant³⁵. The phonon frequency for each K-point is plotted in Fig. 5. As can be seen in Fig. 6, the C_v of NaYRuH₆ and Mn-based AA${}^{\prime}$MH₆ materials are plotted with broken lines, and that of the rest of the materials is indicated by solid lines. It is clear that the C_v of NaYRuH₆ and Mn-based AA${}^{\prime}$MH₆ materials are smaller than that of the other materials at the low temperatures (0 to 150 K). At the low-temperature limit (T ≥ θ_D, θ_D/T < < 1), C_v is proportional to (T/θ_D)³. Since θ_D is higher compared to that of other materials, C_v is smaller at low temperatures for NaYRuH₆ and Mn-based AA${}^{\prime}$MH₆.

Fig. 6: The specific heat capacity (C_v) of AA${}^{\prime}$MH₆ materials as a function of temperature (T).

Here, broken lines indicate C_v of AA’MnH₆ and NaYRuH₆ materials, while solid lines indicate that of the other materials.

Electronic Properties

As can be seen in Table 6, A, A${}^{\prime}$ and M elements lose electrons (except in Ru, where it has small negative value), while H atoms gain electrons. Thus, we can expect an ionic character in A-H, A${}^{\prime}$-H, and M-H bonds. Even though A and A${}^{\prime}$ sites are symmetrically equivalent, the atoms at those sites can lose a different amount of electrons. This is mainly because atoms at those sites have different oxidation states. Based on Table 6, Na, K, and Cs alkali atoms have their usual oxidation state (+1), while alkaline earth atoms such as Ca, Sr, and Ba lose more than 1 electron as they can donate up to 2 electrons. Al, Si, and Tl exhibit their most common oxidation states, which are +3, +4, and +1, respectively. It is reported that first-principles computations provide only negligible changes in the local transition-metal charge for semiconducting crystals³⁶. Therefore, we propose that we can consider MH${}_{6}^{n-}$ complex as a single unit since the M-H bond lengths are very short compared to other H-related bonds. n can be found by computing Δq_M + 6 × Δq_H, which is greater than 2 for all the M atoms except for Ni and Si. For those two atoms, n ≈ 1.6. Therefore, we can expect MH${}_{6}^{2-}$ for Si and Ni complexes, while MH${}_{6}^{3-}$ for the rest of the complexes. Kadir et al. suggest that IrH${}_{6}^{3-}$ complexes exist in AA${}^{\prime}$IrH₆ semiconductors²⁷. Therefore, MH${}_{6}^{3-}$ can be the common complex that exists in AA${}^{\prime}$MH₆ materials.

Table 6 The Bader charge transfer (Δq) in electrons for each element of AA${}^{\prime}$MH₆ Materials.

Figures 7 and 8 show the band structures and partial density of states (PDOS) of the AA${}^{\prime}$MH₆ materials. It is clear that all six AA${}^{\prime}$MH₆ materials are semiconductors. The bandgap for each material is mentioned in Table 7. The DFT calculations with PBE exchange-correlation functional underestimate the band gaps due to self-interaction error. It has been shown that the Heyd-Scuseria-Ernzerhof (HSE) screened Coulomb hybrid functional calculations provide reasonable estimation for the band gaps of semiconductors^37,38. HSE06 uses $\frac{1}{4}$ of exact exchange and $\frac{3}{4}$ of PBE exchange. Based on our HSE06 computations, all the AA${}^{\prime}$MH₆ semiconductors can have bandgaps greater than 2.00 eV (see Supplementary Information). The bandgap range of wide-bandgap semiconductors is considered as the range above 2 eV²³. Thus, we can identify that those materials are wide-bandgap semiconductors. As reported by Kadir et al., NaCaIrH₆, NaSrIrH₆, NaBaIrH₆, KSrIrH₆ and KBaIrH₆ have bandgaps between 2.91 and 3.33 eV²⁷ (see Supplementary Information). Wide-bandgap semiconductors are vital for manufacturing optical devices emitting green, red, and UV frequencies and also power devices functioning at higher temperatures^23,24.

Fig. 7: Band structures of AA${}^{\prime}$MH₆ materials.

Fermi energy marks zero energy.

Fig. 8: Partial density of states (PDOS) of AA${}^{\prime}$MH₆ materials.

Fermi energy marks zero energy.

Table 7 The band gap, conduction band minimum (CBM), valence band maximum (VBM) and type of semiconductor for AA${}^{\prime}$MH₆ materials.

Other than in BaCsAlH₆ and CsKSiH₆, all the AA${}^{\prime}$MH₆ materials have their conduction band minimum (CBM) at X high-symmetric K-point. The CBM of BaCsAlH₆ and CsKSiH₆ are at Γ points. The valence band maximum (VBM) of BaNaRhH₆, SrTlIrH₆, YMnZnH₆, NaYRuH₆, and AgZrMnH₆ exist at W K-point. BaSrZnH₆, KNaNiH₆ and BaSrZnH₆ have VBM at X, while that of CaScMnH₆ and AgZrMnH6 is at K high-symmetric point in the reciprocal space. Thus, both CBM and VBM of BaSrZnH₆ and KNaNiH₆ reside at X K-point, indicating those materials are direct bandgap semiconductors. Direct bandgap semiconductors are preferred for LED and laser devices over their indirect counterparts. Wide-bandgap semiconductors with direct bandgap are widely investigated for solar cells due to optical transparency³⁹. BaNaRhH₆, KNaNiH₆, CaCsMnH₆, and NaYRuH₆ materials have very flat bands near the Fermi level, which is indicated by zero energy. Relative to other materials, BaSrZnH₆ contains narrow (less flat) bands near the Fermi level. As a result, this can lower the effective mass of the carriers. Some research has shown that low effective mass will help developing efficient thermoelectric devices^40,41,42. As shown by electronic band theory, the electron effective mass can be very high in the flat bands⁴³. It is also shown that flat bands at the bottom of the conduction bands can provide high thermoelectric power⁴⁴. YMnZnH₆ and ScZnMnH₆ materials also exhibit that the CBM are relatively flat. Moreover, as shown in Fig. 7, we can modulate the shape of the bands near the Fermi level using the chemical formula. As a result, the thermoelectric properties can be tuned. Therefore, we propose that AA${}^{\prime}$MH₆ semiconductors should be investigated for thermoelectric applications. Our partial density of states (PDOS) studies reveal that d-orbitals of transition metal atoms reside at the M site dominate in the valence region near the Fermi level. Even though the transition metal atoms can be found at A and A${}^{\prime}$ sites, their pdos of d-orbitals are not significant near the Fermi level.

Method

Generative adversarial network

The hypothetical materials used in our research are generated by our CubicGAN²², a generative adversarial network (GAN) based model for generating cubic crystal structures in a high-throughput manner. Our GAN model consists of a generator network and a discriminator/critic network. The discriminator learns to tell real materials from fake materials generated by the generator. The generator learns how to generate samples with similar distribution as the training samples. After trained, we can sample from the generator to generate nonexisting materials. In CubicGAN, we focused on generating ternary and quaternary materials with the space groups 221, 225, and 216. Moreover, to simplify the problem, CubicGAN uses special fractional coordinates, all in the set of {0.0, 0.25, 0.5, 0.75}. The CubicGAN is trained using material data from OQMD^45,46 and is evaluated on material data from Materials Project⁴⁷ and ICSD⁴⁸. The main framework of CubicGAN and the post-processing for the generated materials are shown in Fig. 9. It is notoriously hard to train the original GAN model because the adversarial loss is not continuous in the generator, which causes vanishing gradients and saturation in the discriminator. We take advantage of the Wasserstein GAN with gradient penalty by penalizing the norm of gradients of the critic with respect to the inputs¹⁸. The critic takes real materials and fake materials generated by the generator and then outputs a score which can be interpreted as how real the input materials are. The score is used to update the parameters of the models of the generator and the critic. The adversarial loss is defined as:

$${{{\mathcal{L}}}}=\mathop{\mathbb{E}}\limits_{\tilde{{{{\bf{x}}}}}\sim {{\mathbb{P}}}_{{{{\rm{g}}}}}}[D(\tilde{{{{\bf{x}}}}})]-\mathop{\mathbb{E}}\limits_{{{{\bf{x}}}}\sim {{\mathbb{P}}}_{{{{\rm{r}}}}}}[D({{{\bf{x}}}})]+\lambda \mathop{E}\limits_{\hat{{{{\bf{x}}}}}\sim {{\mathbb{P}}}_{\hat{{{{\bf{x}}}}}}}[{({\parallel {\nabla }_{\hat{{{{\bf{x}}}}}}D(\hat{{{{\bf{x}}}}})\parallel }_{2}-1)}^{2}]$$

(6)

where D means the score function from the critic. $\hat{{{{\bf{x}}}}}$ is the linear interpolation between a real material x and the generated one $\hat{{{{\bf{x}}}}}$ and $\mathop{E}\limits_{\hat{{{{\bf{x}}}}}\sim {{\mathbb{P}}}_{\hat{{{{\bf{x}}}}}}}[{({\parallel {\nabla }_{\hat{{{{\bf{x}}}}}}D(\hat{{{{\bf{x}}}}})\parallel }_{2}-1)}^{2}]$ is the gradient penalty which enforces gradients with the norm at most 1 everywhere. λ is set 10 by default in this work.

Fig. 9: The main framework of CubicGAN.

a WGAN architecture and b postprocessing of generated samples. Z is the random noise. E is the element vector embedded by 23 element properties. SGP is the one hot embedding for the space group.

Conditioning on random noise, three or four-element combinations, and space group, the generator not only generates materials with existing prototypes but also generates stable ones with nonexisting prototypes. When the CubicGAN generates 10 million materials, it can rediscover most of the cubic materials in Materials Project and ICSD. In CubicGAN, we only focus on the generated materials with prototypes, which are defined by the anonymous formula and the space group ID. In total, 24 and 1 nonexisting prototypes are found in 10 million generated ternary and quaternary materials, respectively. Sub-figure (a) of Fig. 9 shows how to filter out the materials. On average, 90% of generated materials have readable CIFs, and we only select materials with neutral charge and negative formation energy predicted by CGCNN⁴⁹. After filtering down materials with nonexisting prototypes, we performed DFT calculations, and 36847 candidate materials have been relaxed successfully. Further, 506 stable materials are verified by phonon dispersion.

Nonmetal - metal classifier

To develop a nonmetal - metal classifier, we first collected the pretty formulas, Bravais lattice type, and bandgap details of all the cubic quaternary materials from the MP database. There were 2578 nonzero bandgap materials (semiconductors and insulators) and 1,438 metals in the collected dataset. We considered 55 elemental and electronic structure attributes, such as the first ionization energy, atomic volume, electronegativity, total number of valence electrons, and number of valence electrons in s, p, d, and f orbitals, to develop the feature set (see Supplementary Information). The weighted average (Avg.) and a maximum difference of those properties for a given chemical formula were added to the feature set. The Avg. of a property S of a quaternary compound A_αB_βC_γD_δ was calculated based on the following expression,

$${S}_{{{{\mbox{A}}}}_{\alpha }{{{\mbox{B}}}}_{\beta }{{{\mbox{C}}}}_{\gamma }{{{\mbox{D}}}}_{\delta }}^{{{{\rm{Avg}}}}}=\frac{1}{\alpha +\beta +\gamma +\delta }(\alpha {S}_{{{{\rm{A}}}}}+\beta {S}_{{{{\rm{B}}}}}+\gamma {S}_{{{{\rm{C}}}}}+\delta {S}_{{{{\rm{D}}}}}),$$

(7)

where S_A, S_B, S_C and S_D are the property S of A, B, C, and D elements, respectively. Altogether, 119 features were considered for training the models.

We created the DNN classifier with two hidden layers using Keras⁵⁰ on top of TensorFlow⁵¹. The first and second hidden layers of DNN include 200, and 100 neurons, respectively. To include the nonlinearity in the system, we shifted the summed weighted inputs of each layer through the rectified linear unit (ReLu) activation function. We randomly dropped out 5% of the units of the hidden layers while training the models. This process is very important for limiting the overfitting of training data. Another useful approach to diminishing overfitting is weight regularization. We employed Ridge (L2) regularization method for adding penalties during updating weights. The adaptive moment estimation (Adam) optimizer with a 0.001 learning rate was considered with binary cross-entropy as the loss function and the metric during the calculations. The optimized number of epochs and batch size are 500 and 1500, respectively.

We developed a random forest classifier (RFC) as the second model, which uses an ensemble technique. Here, data is divided randomly, which is known as bagging and carries out training with multiple decision trees. The final prediction is given by averaging the output of all the decision trees. The hyperparameter optimization was performed using GridSearchCV algorithm as implemented in the Scikit-learn code⁵². The optimized number of decision trees, minimum samples split, minimum samples leaf, and maximum depth are 500, 10, 3, and 90, respectively. Furthermore, we used the RFC model to study the feature importance for whole quaternary materials data set. It will help discovering semiconductors in the future.

For both DNN and RFC models, the cubic quaternary materials dataset with 4016 materials was split randomly into 98 % and 2 % as the training and testing sets, respectively. The 10-fold cross-validation with accuracy as the scoring method was performed on the training set. Here, the training set was partitioned into 10 subsets, where 9 subsets were for training the model and the remaining subset was for validating.

Density functional theory (DFT) calculations

Density functional theory calculations were performed as implemented in the Vienna ab simulation package (VASP) code^53,54,55,56. The electron wave functions were described using the PAW pseudopotentials^57,58. The exchange-correlation interactions were treated based on the generalized gradient approximation (GGA) within the Perdew-Burke-Ernzerhof (PBE) formulation^59,60. The energy threshold value of the plane-wave basis was set as 500 eV. In addition, the energy convergence criteria were set to 10⁻⁸ eV, and the force convergence criterion for the ionic steps is set to 10⁻² eV/Å. The Brillouin zone integrations were performed using a dense k-point mesh within the Monkhorst-Pack scheme for the structure optimizations, band structure, density of states, mechanical properties, and phonon calculations. For instance, a 14 × 14 × 14 K-mesh was used for BaNaRhH₆ with 5.5105 Å lattice constant. The 2 × 2 × 2 supercells were employed for obtaining Phonon dispersions using the Phonopy code³⁵. The elastic constants were calculated by employing density functional perturbation theory (DFPT) as implemented in VASP⁶¹. VASPKIT code³¹ was used to obtain the bulk modulus (K), Shear modulus (G), Young’s modulus (Y), and Poisson’s ratio (μ) of the materials based on the Hill method⁶².