Stability diagrams, defect tolerance, and absorption coefficients of hybrid halide semiconductors: High-throughput first-principles characterization

Organic–inorganic hybrid halide perovskites have emerged as one important class of semiconductor materials for various types of optoelectronic applications, such as solar cells and light-emitting diodes, because of their low-temperature solution-based fabrication techniques and excellent optoelectronic properties.^1–6 As a proof of example, the power conversion efficiency (PCE) of the perovskite-based solar cells at lab-scale testing has increased rapidly from 3.8% to 25.2% in the past few years.⁷ This is mainly attributed to their exceptional optoelectronic properties such as tunable bandgaps, high absorption coefficient, long carrier diffusion length and lifetime, strong defect tolerance, and high carrier mobility.^8–14 

Despite the exceptional properties, it is still of urgent demand to search for stable and non-toxic alternative hybrid materials because of low stability and the presence of toxic lead in the halide perovskite.⁶ Some great research efforts have been made to search for new halide perovskite materials using the high-throughput material design approach,^15–24 although nearly all of them are focusing on either purely inorganic perovskites or a limited number of hybrid perovskite structures (cubic structure or its distorted derivatives). This is mainly due to great computational challenges caused by the structural complexity of the hybrid structures in the large-scale first-principles calculations. It is also worth noting that quantitative physical properties of hybrid perovskites cannot be directly derived from their inorganic bulk analogs.

In our recent work, instead of focusing on the prototype perovskite structure, we have carried out a systematic investigation of all the possible hybrid halide materials on the basis of the 24 prototype structures and five relatively small organic cations using the high-throughput computational material design approach.²⁵ By employing a group of combinatorial material descriptors that cover the critical features of electronic band structures, several energetic parameters, structural integrity at room-temperature, and exciton binding energies, we have successfully identified 13 candidates for solar energy conversion and 23 candidates for light-emitting diodes (a total number of 29 hybrid compounds) in five different types of crystal structures out of a quantum material repository of 4507 hybrid halide materials.²⁵ 

In this work, to provide a more detailed guide for the experimental synthesis and characterization, we systemically studied the stability diagrams, defect tolerance, and absorption coefficients for the identified 29 promising hybrid halide materials. To do this, we first calculated thermodynamically stable ranges for all these compounds and considered non-equivalent lattice sites in each compound for building defect structures. A total number of about 5000 defect structures (including neutral and charged) were computed to analyze the defect transition energy levels using high-throughput first-principles calculations. The absorption coefficients were calculated from the dielectric function. For convenience, two representative compounds, MA₂SnI₄ and MA₃Sb₂I₉, were selected for discussing the computational and analysis process. The calculated stability diagrams, defect formation energies and transition levels, and absorption coefficients for all the 29 compounds are presented in the supplementary material.

The automatic framework AFLOW²⁶ based on the Vienna Ab initio Simulation Package (VASP)²⁷ was used for high-throughput first-principles Density Functional Theory (DFT) calculations. The Projector Augmented Wave (PAW) pseudopotentials were used for describing electron–ion interactions,²⁸ and the Generalized Gradient Approximation (GGA) parameterized by Perdew, Burke, and Ernzerhof (PBE) was used for treating electron–electron exchange-correlation functional.²⁹ The van der Waals (vdW) functional DFT-D3 is incorporated to properly describe the long-range dispersion interactions between the organic molecules.³⁰ A k-point grid of 0.05 Å⁻¹ was automatically set for structural relaxation with a convergence tolerance of 0.01 meV/atom, and a denser grid of 0.04 Å⁻¹ was used in static calculations. Other computational settings such as the cutoff energy (the highest value of the pseudopotentials) were managed by the AFLOW code that generates appropriate entries for the structural relaxation and static calculations sequentially and automatically.²⁶ The supercell model (with more than 40 atoms) was built for each material system to study defect properties, in which the supercell size was determined based on the convergence tests for the formation energy of iodine vacancy; see Fig. A1 and Table A1 in the supplementary material. Our calculations indicate that the supercell with more than 40 atoms is large enough to produce converged defect formation energy. Similar conclusions can also be found in the previous literature on perovskite oxides.^31,32 The hybrid DFT calculations within the Heyd–Scuseria–Ernzerhof (HSE) formalism with 25% Hartree–Fock (HF) exchange are employed to calculate the absorption coefficients.^33,34

All the possible intrinsic point defects, including vacancies (V_A, V_B, and V_X) and interstitial (A_i, B_i, and X_i) and antisite (A_B, B_A, A_X, X_A, B_X, and X_B) defects, were considered for all the 29 hybrid halide compounds A_xB_yX_z; see the list of compounds in Table I. The defect structures were built based on all the possible non-equivalent lattice sites in each compound. The Python Charge Defects Toolkit (PyCDT) based on the Interstitial Finding Tool (InFiT) was used to determine possible interstitial sites.^35,36 The building process yields a total number of 2160 defect structures. After first-principles total energy calculations, we selected the defect structure with the lowest total energy for the same type of defects for the further investigation of defect transition levels. This process leads to a total number of 348 selected structures, which are the lowest-energy site representatives for 12 intrinsic point defects in each of the 29 compounds and used for charged defect calculations.

For charged defect calculations, we considered defect charge ranges derived from all the possible oxidation states of involved element(s) and the defect type based on the selected 348 defect structures. In terms of defect type, vacancies of atom Aⁿ (n is the formal charge of A in the compound A_xB_yX_z) can have charge states in the range of [−n, +n]; interstitial A^m (m represents all the possible oxidation states of A) can have charge states in the range of [min(m), max(m)]; and antisite A_B have charge ranges of [min(p), max(p)], where p represents all possible oxidation states of A and B.³⁵ According to the above charge ranges, a total number of about 2400 charged defect calculations were to be carried out.

It is noted that, in the charged defects calculations under the periodic boundary conditions, the Coulombic interaction between a charged defect and its periodic images is the dominant source of error for defect energetics. The interaction is inversely proportional to the supercell periodic length. The supercells used in the high-throughput calculations are finitely sized, so the charged defect images cannot be effectively isolated. Thus, the Coulombic interaction cannot be completely eliminated either. To correct supercell finite-size effects in the charged defects calculations, we employed the image charge correction proposed by Freysoldt et al. in our calculations by calling the pymatgen package.^36,37

We first computed stability diagrams of the 29 hybrid halide compounds by determining their thermodynamically stable chemical ranges with respect to their decomposers. Note that the dynamic stability, as a more realistic descriptor for screening novel materials than the thermodynamic stability, has been investigated in our prior study using ab initio molecular dynamics (AIMD) simulations.²⁵ Here, the calculated thermodynamic stability diagrams have two major purposes: (i) to determine the chemical potential ranges for synthesizing the hybrid materials and (ii) to choose appropriate chemical potentials for subsequent defect formation energy calculations. The first representative compound, MA₂SnI₄, a tetragonal compound consisting of layers of corner-sharing SnI₆, was taken as an example to discuss the computational process; see Fig. 1(a). Under thermodynamic equilibrium growth conditions, the formation of MA₂SnI₄ should satisfy the following equation:

in which Δμ_X = μ_X − E_X is the chemical potential change in the X component (organic molecule MA was treated as one component) and ΔH is the formation enthalpy. The total energy of MA species (E_MA) was calculated using a standalone neutral MA molecule in a big box, which represents the gaseous phase as the reference ground state for the species. In addition, to prevent the formation of secondary phases MAI, SnI₂, and SnI₄, the following stability limits should also be satisfied:

By applying these thermodynamic conditions, we plotted the stability diagram for MA₂SnI₄ against Δμ_Sn and Δμ_I, as shown in Fig. 1(b). The green region is the thermodynamically stable range for MA₂SnI₄. In this region, we have selected three representative points: A (Δμ_Sn = −4.41 eV, Δμ_I = 0 eV), B (Δμ_Sn = −1.56 eV, Δμ_I = −0.75 eV), and C (Δμ_Sn = 0 eV, Δμ_I = −2.2 eV), which represent chemical potential conditions of X-rich/B-poor, X-moderate/B-moderate, and X-poor/B-rich, respectively.

Similarly, the stability limits and stability diagram of the second example compound MA₃Sb₂I₉ were calculated and plotted, as shown in Figs. 2(a) and 2(b). The three representative Δμ points for MA₃Sb₂I₉ are A (Δμ_Sb = −2.73 eV, Δμ_I = 0 eV), B (Δμ_Sb = −1.04 eV, Δμ_I = −0.35 eV), and C (Δμ_Sb = 0 eV, Δμ_I = −0.91 eV). The calculated stability diagrams for all the 29 compounds, plotted against the chemical potential change Δμ_B and Δμ_X, are shown in Figures B1–B5 of the supplementary material. To deduce all the possible thermodynamic stability limits for each of the 29 compounds, we have carefully considered all the possible neighboring secondary phases on the basis of the open quantum material repositories AFLOWLIB³⁸ and Materials Project.³⁹ The binary halide compounds BX_n for B elements in the 29 A_xB_yX_z compounds were systemically recalculated in the DFT-D3 approach.

Defect tolerance significantly influences the optoelectronic properties of semiconductor materials and can be evaluated from the defect transition levels and the defect formation energy. A strong defect tolerance can prevent the undesired recombination of photoexcited electrons and holes, and it has the following features: the intrinsic defects with a low formation energy will not create deep gap states, while those producing deep gap states have a relatively high formation energy.^9,25 Formation energy (E^f) of a defect X at the charge state q can be calculated using the following equation:⁴⁰ 

in which E_tot[X^q] and E_tot[bulk] are total energies of defected and pristine supercells, respectively. ∑n_iμ_i is the sum of chemical potentials to compensate in the defected supercell, E_F is the Fermi energy, and E_corr is a correction term for electrostatic interactions between supercell images. Here, we only considered defect charge states derived from involved elements’ oxidation states in the compound, which gives upper and lower limits of defects’ charge states. All the charge states between these limits are plotted in the formation energy plots as a function of Fermi energy. This consideration leads to a convenient visualization of the most probable charge states and avoidance of redundancy in the plots, which is also well justified by a prior computational study of point defects in the hybrid halide systems.⁹ By taking iodine as one example, one only needs to consider its oxidation state of 1− in the hybrid iodides. As a result, an interstitial iodine (I_i) can have two possible charge states of 0 and 1-, and one iodine vacancy (V_I) can have two possible charge states of 0 and 1+. The reason why iodine’s positive oxidation states are not considered is that there are no elements oxidizing iodine in the hybrid metal iodide materials.

Figures 1(c) and 2(c) show the formation energies of intrinsic point defects under chemical potential conditions A, B, and C for MA₂SnI₄ and MA₃Sb₂I₉, respectively. The lower E^f are shown in colored solid lines, while higher E^f are shown in black dashed lines. The slope of the line represents charge q, so horizontal lines show E^f of neutral defects (q = 0). E^f lines change as the ∑n_iμ_i term changes among A, B, and C, which is in agreement with chemical condition changes in the involved elements. For example, as the chemical condition becomes more I-rich, E^f of V_I and I_i increases and decreases, respectively. Accordingly, the intersection points between different charge states of a defect correspond to the same E_F. As shown in Fig. 1(c), the defects with relatively low E^f in MA₂SnI₄ are V_Sn, I_Sn, and MA_Sn at A; V_I and Sn_i at B; Sn_I, MA_I, MA_i, V_I, and Sn_i at C.

The transition energy level (ε(q1/q2)) of a defect is the Fermi level position that corresponds to the intersection point of two different charge states (q1 and q2) of a defect in the E^f plots, where the formation energies of the two charge states are equal. In other words, ε(q1/q2) is the energy level where the defect can accept or donate electrons so that the charge states of a defect can be changed. As shown in Fig. 1(c), the lowest E^f of V_Sn (green line) is composed of two parts, i.e., $VSn0$ (left) and $VSn2−$ (right), and the two parts intersect at E^f = 1.02 eV. At this transition energy level, $VSn0$ accepts two electrons and becomes $VSn2−$⁠. On the basis of Eq. (5), ε(q1/q2) can be calculated using the following equation:⁴⁰ 

where E^f(X^q; E_F = 0) is the formation energy of defect X^q, when the Fermi level is at 0 [valence band (VB) maximum]. One defect with ε(q1/q2) close to the middle of the bandgap has a deep transition level that will attract electrons/holes and acts as undesired nonradiative recombination centers. In contrast, defects with ε(q1/q2) close to band edges or at the inside of the conduction band (CB) or valence band (VB) will not significantly affect the optoelectronic performance.

The calculated transition levels of MA₂SnI₄ and MA₃Sb₂I₉ are shown in Figs. 1(d) and 2(d), respectively. In the ε(q1/q2) figures, the red and blue levels indicate intrinsic donors and acceptors, and the top and bottom gray regions indicate CB and VB, respectively. For MA₂SnI₄ (a bandgap of 1.42 eV), as shown in Fig. 1(d), the defects that induce deep transition levels are Sn_I (ε(1+/2+) = 1.04 eV and ε(2+/3+) = 0.52 eV), MA_I (ε(0/1+) = 1.03 eV and ε(1+/2+) = 0.22 eV), V_Sn (ε(0/2−) = 1.02 eV), I_MA (ε(0/1− = 1.06 eV), and V_MA (ε(0/1−) = 0.74 eV). Among these defects, I_MA has very high formation energies at all chemical conditions, as indicated by the dashed lines, and, thus, will not lead to deep gap states, indicating a high defect tolerance. Sn_I (yellow) has a high formation energy at chemical conditions A and B, MA_I (orange) has a high formation energy at chemical conditions A and B, V_Sn (green) has a high formation energy at chemical conditions B and C, and V_MA (lime) has a relatively high formation energy at all chemical conditions; see Fig. 1(c). Therefore, the best chemical condition to synthesize MA₂SnI₄ is at B (I-moderate/Sn-moderate), under which undesired Shockley–Read–Hall (SRH) nonradiative recombination centers can be largely prevented in the material.

For MA₃Sb₂I₉ (bandgap= 2.02 eV), as shown in Fig. 2(d), the defects that create deep transition levels are Sb_I [ε(2+/4+) = 0.86 eV] and Sb_i [ε(1+/2+) = 1.10 eV and ε(2+/3+) = 0.24 eV]. Both Sb_I (yellow) and Sb_i (brown) have high formation energies at chemical condition A, but very low formation energies at B and C. Therefore, to prevent the deep transition levels in MA₃Sb₂I₉, the material should be synthesized under chemical condition A (I-rich/Sb-poor).

The calculated defect formation energies and transition energy levels for all the 29 compounds are shown in Figs. C1–C29 of the supplementary material. To determine the defect tolerance for each compound, we listed the defects that have low formation energies and, meanwhile, can create deep transition levels in Tables D1–D5 of the supplementary material. The defect tolerance was determined based on the three chemical conditions A, B, and C and also summarized in Table I. In specific, for one compound, if there is at least one chemical condition that prevents all kinds of defects mentioned above, the compound is determined to have high defect tolerance; if such defects cannot be prevented at any chemical condition, the compound is determined to have low defect tolerance. From this perspective, there are 15 out of 29 compounds that show high defect tolerance.

In this section, we investigated the optical absorption coefficients using the following equation:⁴¹ 

where ω is the photon frequency and c is the speed of light in vacuum. The real part ε₁(ω) and imaginary part ε₂(ω) are from the complex dielectric function ε(ω). The real part ε₁(ω) follows the Kramer–Kronig relationship, and the imaginary part ε₂(ω) is calculated from the appropriate momentum matrix elements between the occupied and the unoccupied wave functions within the selection rules over the Brillouin zone. On the basis of the ground-state electronic structure calculations at the HSE06 level, we obtained ε₂(ω) and ε₁(ω) of the dielectric function and calculated the optical absorption coefficient of the screened 29 compounds.

The calculated absorption coefficients as a function of wavelength (nm) for the two representative compounds, MA₂SnI₄ and MA₃Sb₂I₉, are shown in Figs. 3(a) and 3(b), respectively. The averaged (non-polarized) absorption coefficients at 450 nm are estimated to be 2.4 × 10⁵ cm⁻¹ for MA₂SnI₄ and 3.1 × 10⁵ cm⁻¹ for MA₃Sb₂I₉, even higher than 2 × 10⁵ cm⁻¹ for MAPbI₃. Interestingly, MA₂SnI₄ shows a much higher absorption coefficient along the in-plane (x and y) direction than that along the out-of-plane (z) direction; see Fig. 3(a). The anisotropy of the calculated absorption coefficients resembles that of the two-dimensional Ruddlesden–Popper perovskites such as Ge-based compound in our prior work, which is attributed to the two-dimensional structural feature.¹⁴ However, it is noted that in these quantum well-like perovskites, local fields and excitonic effects have major impacts on their optical properties. To accurately describe optical properties of these materials, DFT calculations with GW approximation and the Bethe–Salpeter equation (BSE) are required,^42,43 which, however, is more time-consuming. Therefore, to have a direct comparison with the well-studied MAPbI₃ at the same computational level, the ground-state HSE06 calculations were adopted in this work. All the calculated absorption coefficients for the 29 compounds are shown in Figs. E1–E5 of the supplementary material.

In summary, by using high-throughput first-principles electronic structure calculations, we have systematically studied the stability diagrams, defect tolerance, and absorption coefficients for the screened 29 lead-free candidate hybrid semiconductors for optoelectronic applications. The calculated stability diagram outlines the thermodynamically stable range for the equilibrium growth of the predicted compound with different chemical potentials. The defect tolerance was evaluated from the calculated defect formation energies and transition levels, and the absorption coefficients were calculated from dielectric functions. These computational studies provide a detailed guide to the further experimental synthesis and characterization of these hybrid compounds, with a potential to facilitate the development of novel optoelectronic devices.

See the supplementary material for the complete phase diagram, defect formation energy, transition energy level, defect tolerance, and absorption coefficient for all 29 candidate compounds studied.

This work was supported by the Academic Senate General Campus Research Grant Committee at the University of California, San Diego, the National Science Foundation under Award No. ACI-1550404, and the Global Research Outreach (GRO) Program of Samsung Advanced Institute of Technology under Award No. 20164974. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation, Grant No. ACI-1548562. D.M.-L. and V.R.-V. acknowledge the ENLACE bi-national summer research program at the University of California, San Diego.