Computational Screening of All Stoichiometric Inorganic Materials.

Forming a four-component compound from the first 103 elements of the periodic table results in more than 1012 combinations. Such a materials space is intractable to high-throughput experiment or first-principle computation. We introduce a framework to address this problem and quantify how many materials can exist. We apply principles of valency and electronegativity to filter chemically implausible compositions, which reduces the inorganic quaternary space to 1010 combinations. We demonstrate that estimates of band gaps and absolute electron energies can be made simply on the basis of the chemical composition and apply this to the search for new semiconducting materials to support the photoelectrochemical splitting of water. We show the applicability to predicting crystal structure by analogy with known compounds, including exploration of the phase space for ternary combinations that form a perovskite lattice. Computer screening reproduces known perovskite materials and predicts the feasibility of thousands more. Given the simplicity of the approach, large-scale searches can be performed on a single workstation.

Introduction

Currently, over 184,000 entries in the Inorganic Crystal Structure Database (ICSD) involve 9,141 structure types; 66,814 of these materials have also been subject to quantum mechanical calculations, and information on their basic electronic structures and thermodynamics is included in the Materials Project (powered by the PYMATGEN infrastructure).

The configurational phase space for new materials is immense, but blind exploration of the periodic table is a daunting task. Fortunately, over a century of research in the physical sciences has provided us with myriad rules for assessing the feasibility of a given stoichiometry and the likelihood of particular crystal arrangements. Examples of chemical phenomenology include the radius ratio rules and Pettifor maps for structure prediction, as well as electronegativity and chemical hardness for predicting reactivity. Pauling's rules provide predictive power for ionic or heteropolar crystals. A wealth of knowledge exists for understanding the physical properties of tetrahedral semiconductors. Recent examples of searches for new materials that draw from existing chemical knowledge include 18-electron ABX compounds, hyperferroelectric superlattices, and organic-inorganic perovskites.11, 12

The reliability and predictive power of atomistic material simulations is increasing.13, 14 Many approximations are being removed as high-performance supercomputers reach petaflop scale. This includes more accurate quantum mechanical treatment of electron-electron interactions in the solid state, as well as more realistic models of chemical disorder. However, because of the computational cost, high-throughput screening with first-principle techniques is usually limited to hundreds or thousands of materials—a small fraction of the overall phase space.

We report a procedure for navigating the materials landscape with low computational effort, and it can be achieved with simple chemical descriptors. We first explore the magnitude of the task at hand by enumerating combinations of elements and ions for binary, ternary, and quaternary compositions. We demonstrate that chemical constraints can narrow the search space drastically. Examples of how deeper insights can be gained are illustrated for electronic (photoelectrodes for water splitting) and structured (perovskite-type) materials. The procedure can be used to comfortably explore the vast compositional space or as the first step in a multi-stage high-throughput screening process. Instead of being a roadblock to achieving new functionality, the combinatorial explosion for multi-component compounds provides fertile ground for the discovery of innovative materials.

Results

Elemental Combinations

To begin, one can map chemical space by enumerating the ways in which the constituent elements of the periodic table can combine. If we restrict ourselves to the first 103 elements (to the end of the actinide series), the combinations (i.e., ) for two, three, and four components are 5,253, 176,851, and 4,421,275, respectively. For five components, the combinations exceed 87 million.

Physically, the situation is more complex. Elements can combine in different ratios, leading to variation in material stoichiometry, e.g., the binary combinations AB, AB2, A2B3, and A3B4. Given elements can also adopt multiple oxidation states, each with a unique chemical behavior, e.g., Sn(II)O, Sn(IV)O2, and Sn(II)Sn(IV)O3. For our enumeration of feasible compounds, we next consider the accessible oxidation states of each element in stoichiometry up to quaternary ABCD, where the integers w, x, y, and z ≤ 8. This definition includes, for example, common ternary pyrochlore oxides (A2B2O7) and quaternary double perovskites (A2BCO6). Using the most common oxidation states extends the first 103 elements of the periodic table to 403 unique ions.

The number of combinations is now drastically increased, as shown in Table 1, such that four-component candidate materials exceed 1012. In order to reduce this composition space, we can introduce selection rules (filters) from chemical theory.

Table 1

Estimates for the Number of Possible Inorganic Materials Allowing for Variable Oxidation States and Stoichiometry with the Constraints of Charge Neutrality and Electronegativity Balance

Type	Constrainta	Number
AwBx	–	3,483,129
AwBx	q	58,614
AwBx	q + χ	14,721
AwBxCy	–	4,753,229,039
AwBxCy	q	174,081,685
AwBxCy	q + χ	32,157,899
AwBxCyDz	–	4,139,315,402,300
AwBxCyDz	q	267,381,955,246
AwBxCyDz	q + χ	32,381,953,858

q, charge neutralitiy; X, electronegativity balance.

We note that the estimations discussed here represent a lower limit on the number of accessible materials. We consider regular inorganic compounds and exclude, for example, non-stoichiometry, organic systems, hybrid organic-inorganic materials, electrides, and intermetallics, for which additional considerations are required for predicting viability.17, 18, 19

Chemical Filters

Rule 1: Charge Neutrality

Ions tend to combine into charge-neutral aggregates. The same thinking applies to both ionic solids and more covalently bonded semiconductors. Any periodic solid must be charge neutral; otherwise, it would have an infinite electrostatic potential. Balancing oxidation states and fulfilling the valence octet rule are equivalent, e.g., III–V semiconductors, such as GaAs, can be represented as Ga3+As3−. Our implementation is inspired by the work of Pamplin and Goodman on the subject of multi-component semiconductors.

A charge-neutrality constraint significantly reduces the total number of candidate materials. The rule states that the formal charges (q) of the components sum to 0, i.e.,

Charge neutrality contracts the compositional space by at least an order of magnitude for binaries, ternaries, and quaternaries (Table 1).

Rule 2: Electronegativity Balance

Further to assuming that all charge-neutral combinations of oxidation states are accessible, we can implement a second constraint based on the electronegativity of the component elements. The empirical electronegativity (χ) scale represents the “attraction” of a particular atom for electrons. For a stable compound, the relation χcation < χanion should be obeyed, i.e., the most electronegative element carries the most negative charge. Here, we employ the Pauling electronegativity scale, which reduces the allowed compositions by a factor of between 4 and 10 for the different numbers of components (Table 1).

It is also instructive to consider existing materials databases (the ICSD and Materials Project). For binary compounds, we find fewer combinations from our estimates as implemented in SMACT (Semiconducting Materials from Analogy and Chemical Theory) than from the ICSD (Figure 1), which can largely be attributed to our exclusion of intermetallics and polymorphs. In the Materials Project, multiple entries for a single crystal structure and chemical composition are removed, and the number of compositions are in close agreement. For ternaries and quaternaries, the compositions passing both charge and electronegativity tests continue to rise exponentially, whereas the number in existing databases remains relatively constant. The increased complexity of ternary and quaternary systems means that their synthesis, characterization, and reporting are more challenging than for binary systems. Nevertheless, the large differences between the numbers of potential and reported materials suggest that wide areas of unexplored compositional space may contain stable and useful materials.

Figure 1

Counting the Number of Possible Multi-component Materials

(Left) Narrowing of compositional space for inorganic materials by imposing chemical constraints of charge (q) and electronegativity (χ).

(Right) Comparison of the accessible materials predicted by SMACT and those reported in the ICSD and the Materials Project.

The numbers reported in this section are vast, and using modern electronic-structure techniques to perform quantitative screening for application is unimaginable. Exploration of the hitherto neglected compositional space will require further guidelines. In the following sections, we demonstrate how additional descriptors can be applied for identifing materials for specific applications.

Compositional Descriptors

Several useful properties can be estimated from knowledge of the chemical composition alone, and here we explore the application of some of these approaches.

Descriptor 1: Electronic Chemical Potential

The concept of atomic electronegativity has been successfully extended to solids, where the geometric mean becomes the single-value descriptor, i.e.,

This descriptor represents a mid-gap energy between the filled (valence band) and empty (conduction band) electronic states. This corresponds to the electronic chemical potential (Fermi level) at 0 K. Butler and Ginley found a linear correlation between the solid electronegativity and the electrochemical flat-band potentials for a range of semiconductors. This was subsequently extended to a wider dataset including metal oxides, chalcogenides, and halides. The method provides a rapid procedure for the estimation of absolute electron energies for multi-component materials. It is now commonly used in the computational screening of new materials for electrochemical applications.25, 26, 27, 28

Descriptor 2: Electronic Structure

Many tight-binding model Hamiltonians exist for semiconductors and dielectrics. One recent approach is based on the atomic solid-state energy (SSE) scale, which provides information on valence and conduction bands on the basis of the frontier orbitals of the constituent ions. Whereas the Mulliken definition of electronegativity is an average of the ionization potential (IP) and electron affinity (EA) of an atom, the SSE reflects the IP of an anion (filled electronic states) and EA of a cation (empty electronic states). The energies of 40 elements were fitted from a test set of 69 closed-shell binary inorganic semiconductors, which has recently been extended to 94 elements. According to the tabulated SSE scale, the band gap (Eg) can be estimated as

For multi-component systems, the limiting values (cation with highest EA and anion with lowest IP) are used. The SSE has a root-mean-square deviation of 0.66 eV against the measured band gaps of 35 ternary semiconductors (see Table S1). This simple method allows for rapid screening of band gaps and absolute band-edge alignment.

Both methods (Equations 2 and 3) have been implemented for arbitrary compositions on the basis of tabulated atomic data in the SMACT package. Because no crystal-structure information is included at this level, the results are qualitative, and the models do not distinguish, for example, between polymorphs.

Electronic Structure: Photoelectrodes

We now use the compositional space and chemical descriptors defined above to search for potential materials for solar fuel generation via photoelectrochemical water splitting.

The properties that are required for viable photoelectrodes include (1) a band gap in the visible range of the electromagnetic spectrum to absorb a significant fraction of sunlight and (2) upper valence and lower conduction bands bridging the water oxidation and reduction potentials, enabling the redox reaction. We set an optimal band-gap range between 1.5 and 2.5 eV. Although the free energy for water dissociation is 1.2 eV, the combination of loss mechanisms found in practical devices could require a band gap as large as 2.2 eV.31, 32

Metal oxides combine many attractive properties for water splitting (e.g., stability and cost), but they usually have band gaps too large to absorb a significant fraction of sunlight. The formation of multi-anion compounds offers a route to modifying the electronic structure. We consider ternary metal chalcohalides (i.e., ABC), with B = [O,S,Se,Te] and C = [F,Cl,Br,I]. We restrict the A cations to those with an SSE higher than the water reduction potential (approximately −4.5 V in relation to the vacuum at pH = 0). The conditions of charge neutrality and electronegativity are used for performing an initial screening that yields 52,094 combinations. With the additional band-gap criterion, the combinations are reduced to 7,676, and the pool of cations is reduced from 25 to 7 with A = [B,Ti,V,Zn,Ga,Cd,Sn]. We further rule out any boron-containing combinations at this stage, because these are known to form discrete molecular units (e.g., BClSe).

Finally, we screen compositions according to the environmental sustainability of the elements. We use the Herfindahl-Hirschman Index (HHIR), which has been developed in the context of thermoelectric applications, for elemental reserves. This index includes factors such as geopolitical influence over materials supply and price. The HHIR for a given composition can be obtained as the weighted average over the constituent elements. At this stage, because stoichiometry is variable, we consider the mean HHIR for each A1B1C1 combination.

The band-edge positions of the 20 candidates with the smallest HHIR values are presented in Figure 2. The HHIR has the effect of eliminating all combinations containing Ga, Te, and Br (although relatively abundant, most of the world's Br is produced from the Dead Sea, making it geopolitically sensitive, as reflected in a high HHIR). There are no entries in the ICSD for most of the candidates that we identified; however, reports can be found for Cd2O6I2, Sn2SI2, and Zn6S5Cl2.34, 35, 36 Both Cd2O6I2 and Sn2SI2 feature in the Materials Project and have band gaps of 3.3 and 1.6 eV, respectively, calculated within density functional theory (DFT). These compare with the SSE band gaps of 2.5 and 2.0 eV. The third compound, Zn6S5Cl2, is reported to have an optical gap of 2.7 eV, which compares with the SSE band gap of 2.4 eV.

Figure 2

Calculated Band-Edge Positions of 20 Promising Element Combinations for Water-Splitting Applications

Band-edge positions were calculated in relation to the vacuum level on the basis of the solid-state energies of the constituent elements. Blue dashed lines indicate the water reduction (above) and oxidation (below) potentials with respect to vacuum.

Only one oxygen-containing compound survived the band-gap screening criterion; the values for metal oxyhalides are generally too large. For OI, the iodide forms the upper valence band (low binding energy of I 5p), whereas it is the oxide (O 2p) for other halides. However, the sensitivity of the oxide ion to its crystal environment is well documented,27, 37 and consequently its SSE carries the greatest uncertainty. This is one aspect where knowledge of the local structure (electrostatic potential) could significantly improve the accuracy of the results.

We must connect composition to crystal structure in order to make more accurate property predictions. Global optimization of crystal structures from first principles is a formidable task, although great progress is being made in this area. We instead adopt an approach based on analogy with known structures through chemical substitutions, as developed by Hautier et al. It uses data-mined probability functions, as implemented in the Materials Project.

To demonstrate the translation from composition to material, we performed crystal-structure mining for the four combinations with the lowest HHIR. The 88 predicted structures were then subjected to a full DFT lattice optimization procedure and ranked by total internal energy. Finally, accurate band gaps were predicted for the lowest-energy structures by hybrid DFT (HSE06 electron exchange and correlation40, 41). The compound Sn5S4Cl2 has an indirect band gap of 1.6 eV and a direct gap of 1.8 eV, which lies within the target range. The band gaps of the other three lowest-energy compounds were calculated to be between 3.0 and 3.4 eV. Full details of the workflow (Figure S1) and band gaps (Table S2) can be found in the Supplemental Information.

The newly identified compound, Sn(II)5S4Cl2, adopts a structure formed from two distinct Sn centered polyhedra: (1) a distorted octahedron with equatorial S and apical Cl ions and (2) a distorted tetrahedron with 4 S ions and a stereochemically active Sn lone pair (Figure S2). The polyhedra form interlocking chains in three dimensions. The electronic density of states reveals an upper valence band composed of hybridized Sn s – Cl p orbitals; such Sn s-based valence bands are considered promising indicators for hole mobility. The lower conduction band is composed mainly of overlapping Sn p orbitals. The chemical structure and bonding characteristics suggest that this material should have favorable carrier transport, crucial for optoelectronic applications.

Crystal Structure: Perovskites

One of the most successful approaches to discovering new materials is structural analogy. The concept is to take a crystal structure with a known chemistry and to replace elements within the structure to tune the physical properties. In the simplest case, this involves direct isovalent substitution, e.g., Zn(II)S → Cd(II)S. Structural analogy can be extended to aliovalent cross-substitution (also termed cation mutation), e.g., Zn(II)S → Cu(I)Ga(III)S2. A systematic methodology was outlined more than 40 years ago in a paper by Pamplin for enumerating charge-neutral tetrahedral semiconductors.

The challenge of going beyond tetrahedral semiconductors is predicting crystal structure. The radius of ions within a lattice has a long history as a geometric descriptor of structural stability. A key example is the application of radius ratio rules by Goldschmidt to predict the propensity of a ternary ABC3 combination to form the perovskite structure:where t is the tolerance factor and r is the ionic radius. Values of t > 1 imply a relatively large A site favoring a hexagonal structure, 0.9 < t < 1 predicts a cubic structure, and 0.7 < t < 0.9 means that the A site is small, preferring an orthorhombic structure. For t < 0.7, other (non-perovskite) structures are predicted. These rules have recently been extended to describe structure-property relationships in hybrid organic-inorganic perovskites.11, 12

In this section, we apply our screening procedure to include knowledge of the crystal structure and estimate the size of the perovskite materials space. We start by enumerating the elemental combinations. We then reduce the set by requiring an octahedral coordination environment for the B site, as contained in the Shannon dataset, and require a combination of oxidation states that are charge neutral. This list is then assessed in terms of t, as defined by the Shannon ionic radii.

We consider single-anion compositions based on C = [O,S,Se,F,Cl,Br,I]. The charge-neutrality and octahedral B-site constraints reduce the 176,851 elemental combinations to 41,725. The tolerance-factor constraint, 0.7 < t < 1.0, further reduces this to 26,567. For potential applications in the energy sector, we can consider candidates with HHIR smaller than that of CdTe (a commercial thin-film photovoltaic material), resulting in a final population of 13,415.

For each anion, an orthorhombic perovskite structure is the most common prediction, and hexagonal is the most rare (Figure 3). The fraction of cubic perovskite structures remains roughly constant within the respective halide and oxide or chalcogenide series; however, it is more dominant for the halides. The presence of Br or I makes a material less sustainable (higher HHIR); otherwise, there is little to differentiate the anions.

Figure 3

Counting Experiments with Perovskites

(A) Combinations found at each stage of the screening procedure.

(B) Perovskite compounds with an HHIR lower than CdTe (3,296) for each anion.

(C) The distribution of hexagonal, cubic, and orthorhombic perovskite structures predicted on the basis of the Goldschmidt tolerance factor and Shannon radii of the ions.

(D) ABC3 combinations found in the Materials Project database. They are sorted into structure type according to space group (here, orthorhombic and lower-symmetry perovskites are grouped together).

Far more oxide and chalcogenide perovskites are predicted than halides. The higher anion charge allows for three distinct cation combinations (I-V, II-IV, and III-III), whereas halides have only I-II. In addition, a greater radius compatibility is found for the group VI anions. We find that the number of plausible perovskite structures increases with the anion radius; however, the lower crustal abundance for heavier elements reduces the number that meet the sustainability criterion.

A search of the Materials Project over the same anion space reveals 920 materials, a small fraction of those predicted from SMACT (26,567). The search includes all standard perovskite space groups. For oxide perovskites, 8.26% of those identified from SMACT are found in the Materials Project; for sulfides, this falls to 0.45% and to 0.12% for selenides. To some extent, the greater number of oxide perovskites discovered reflects the greater research activity in this field; however, synthesis of chalogenide perovskites has been reported,46, 47, 48 and there is interest in these materials for technological applications.49, 50 Of the ABC3 materials reported in the Materials Project, 48% of oxides, 35% of sulfides, and 20% of selenides are in perovskite space groups.

Why are there so few chalcogenide perovskites? The tolerance-factor arguments that work well for metal oxides may not hold for chalcogenide perovskites. Oxygen forms more ionic compounds because of its higher electronegativity and lower polarizability than those of S, Se, and Te. When covalent bonding becomes prevalent, it is known to result in deviations from tolerance-factor behavior. An example is the case of NaSbO3, for which t = 0.92 is commensurate with the formation of cubic perovskite but which forms the non-perovskite ilmenite structure. Goodenough and Kafalas explained this deviation as a result of strong σ bonding between Sb and O.

This procedure demonstrates the power of searching through materials on the basis of structural analogy. Only a small fraction of possible perovskite materials have been synthesized. Although some may not represent thermodynamic ground states, they could be accessible through kinetic control of crystal growth or templated on a substrate. Many interesting chalcogenide perovskites are waiting to be discovered. The final pool of 13,415 feasible compositions is within the grasp of explicit computation by quantum mechanical methods, albeit as part of an ambitious project. Indeed, high-throughput screening of 5,400 multi-anion cubic perovskite structures via DFT has been reported25, 52 and revealed 32 promising new materials for water-splitting applications.

Conclusions

We have demonstrated the utility of chemical theory in quantifying the magnitude of the compositional space for multi-component inorganic materials. Even after the application of chemical filters, the space for four-component materials exceeds 1010 combinations. We further estimate that the five-component space exceeds 1013 combinations. There are many applications in which materials with even higher-order compositions have been developed, e.g., in high-temperature superconductors, where six to seven component materials are common. The number of potential materials is not infinite, but it is immense. The scale of the combinatorial explosion emphasizes the need for effective material-design procedures that employ existing chemical and physical knowledge in a targeted manner. Stochastic sampling of this chemical space is unlikely to be effective in yielding materials with specific functionality. We have presented a procedure that uses simple descriptors to support materials exploration, discovery, and design.

Experimental Procedures

Code Implementation

SMACT, the Python toolkit developed in this work, is available online at https://github.com/WMD-group/SMACT. It is free software made available under the Gnu Public License (GPL) version 3.

All element counts and plots presented in this paper were created with custom codes based on SMACT and written in the Python 2.7 programming language. Elemental data are collated from multiple sources (see Table 2) and made algorithmically accessible in a unified object-orientated interface. Example routines that check element and oxidation-state combinations against the conditions of charge neutrality and electronegativity are provided.

Table 2

Data Sources for SMACT

Data Type	Source
Abundance	estimated crustal abundance of elements from the CRC Handbook of Physics and Chemistry53
Atomic mass	NIST Standard Reference Database 144;54 where the relative abundance of isotopes was unknown or a range of values was provided, a simple mean was taken
Covalent radius	scientific paper55
Electron affinity	scientific paper;56 no default value was used for elements that lack data on electron affinity
Eigenvalues	highest occupied p-state and s-state eigenvalues were tabulated by Harrison57 from the approximate Hartree-Fock calculations of Herman and Skillman58
HHI	elemental Herfindahl-Hirschman Index calculated from geological and geopolitical data33
Ionization potential	NIST Atomic Spectra Database59
Pauling electronegativity	updated values of electronegativity on Pauling's scale were compiled in the CRC Handbook;53 for elements 95 (Am) and above, Pauling's recommended value of 1.3 was used;60 the value for krypton (3.0) was derived from the bond energy of KrF2 and reported in a scientific paper61
SSE	“solid-state energy” model of semiconductors and dielectrics29, 30
SSE (Pauling)	extended estimates of solid-state energy from the correlation between known values and Pauling electronegativity30

Where possible, values recommended by the National Institute of Standards and Technology (NIST) were used.

Scripts that generate the results and plots reported in this paper are available with the SMACT codes. A number of tutorials working through the combinatorial explosion are provided at https://github.com/WMD-group/SMACT_practical.

The codes, collectively named Semiconducting Materials by Analogy and Chemical Theory, are inspired by the pen-and-paper procedure reported by Pamplin in 1964.

Author Contributions

All authors contributed to the development of the SMACT package; the primary coding was performed by K.T.B., A.J.J., and D.W.D. A.W., D.W.D. and K.T.B. wrote the first draft of the manuscript with input, discussion, and analysis from all co-authors.