Predicting the thermodynamic stability of double-perovskite halides from density functional theory

The organic-inorganic hybrid halide perovskite CH₃NH₃PbI₃ (MAPbI₃) has attracted intense attention as a high-efficiency photovoltaic semiconductor. However, the poor stability and toxicity limited its large-scale commercialization, which makes it necessary to search for other new materials as alternatives.^1–6 Among them, double-perovskite halides spurred great attention.^6–17 

The double-perovskite halides are a large family of quaternary halides. The large compositional space of these materials offers opportunities in finding novel optoelectronic materials; nevertheless, it also presents challenges in screening a large number of compounds for identifying the most promising ones with optimal energetic, structural, and electronic properties. Volonakis et al.⁷ designed a series of inorganic double-perovskite halides based on Bi or Sb and noble metals through the density functional theory (DFT) calculations with local density approximation (LDA),¹⁸ and synthesized Cs₂AgBiCl₆ successfully. Filip et al.¹⁵ investigated the electronic structure of Cs₂AgBiBr₆ and Cs₂AgBiCl₆ with LDA. McClure et al.¹⁴ synthesized Cs₂AgBiBr₆ and Cs₂AgBiCl₆, and studied their electronic structure using both Perdew, Burke, Ernzerhof (PBE¹⁹) type of generalized gradient approximation (GGA) and the hybrid functional. Xiao et al.¹⁶ and Savory et al.⁶ studied their thermodynamic stability using PBE and PBE-sol (another type of GGA functional) functionals, respectively. Volonakis et al.⁸ also devised Cs₂AgInX₆ (X = Cl, Br, I) double-perovskite halides and synthesized Cs₂AgInCl₆ successfully, which exhibits a direct bandgap of 2.1 eV according to the hybrid functional calculations based on the crystal structure relaxed by LDA calculations. Interestingly Zhao et al.⁹ explored a series of double-perovskites with a formula A₂BCX₆ (A = Cs⁺; B = Na⁺, K⁺, Rb⁺/Cu⁺, Ag⁺, Au⁺/In⁺, Tl⁺; C = Bi³⁺, Sb³⁺; X = F⁻¹, Cl⁻¹, Br⁻¹, I⁻¹) and predicted 14 stable halides with PBE.¹⁹ Furthermore, Wei et al.¹² reported a hybrid double-perovskite (CH₃NH₃)₂KBiCl₆ and synthesized (CH₃NH₃)₂AgBiBr₆ with an indirect bandgap of 2.0 eV.¹¹ 

Among these theoretical studies, the LDA and GGA (PBE or PBE-sol) functionals were used for relaxing crystal structures and predicting thermodynamic stability. However, it was revealed in previous studies on MAPbI₃ that the van der Waals (vdW) interaction is strong in these perovskite halides and the calculated lattice constants by LDA and GGA functionals deviate significantly from the experimental values.^20–23 When the vdW functionals were employed, the accuracy of theoretically calculated lattice constants could be improved. It was concluded that using vdW functionals is important for the theoretical study of MAPbI₃. For the double-perovskites halides, it is still an open question whether the LDA and GGA functionals can give an accurate description of lattice constants and thermodynamic stability. Therefore, a comprehensive study of the effect of various functionals on the calculated properties is urgently needed for identifying a reliable functional for the high-throughput computational design of new double-perovskites halides.

In this paper, we performed a systematic study on the thermodynamic stability of representative double-perovskite halides Cs₂AgBiCl₆ and Cs₂AgBiBr₆ using six different exchange-correlation functionals, including LDA,¹⁸ two GGA functionals (PBE and PBE-sol),^19,24 and three van der Waals functionals (vdW-DF2, optB86b-vdW, and optB88-vdW).^25–27 We found that the prediction of thermodynamic stability of Cs₂AgBiBr₆ is very sensitive to the choice of the functional. The DFT calculations using the extensively used PBE functional show that Cs₂AgBiBr₆ tends to decompose into binary and ternary compounds and is, therefore, not stable. This result contradicts the successful experimental synthesis of Cs₂AgBiBr₆.¹³ Detailed analysis reveals that the above contradiction results from the failure of PBE calculations to predict the correct ground-state structure of the binary compound AgBr. Moreover, incorrect ground-state structures of a number of binary halides (such as AgCl and CsCl) are predicted by PBE calculations. Extensive tests of a range of different functionals show that LDA, optB86b-vdW, and optB88-vdW functionals can predict the ground-state structures of binary halide compounds such as AgBr, AgCl, and CsCl correctly. Furthermore, optB88-vdW calculations are found to consistently give better results on the calculated lattice constants of all binary, ternary, and quaternary halides compared to other functionals. In view of these two aspects, the optB88-vdW functional is the most reliable functional for predicting the crystal structure and thermodynamic stability of double-perovskite halides. The influence of the functionals on the calculated electronic band structures was also studied. The results show that using different functionals has little influence on the dispersion of the conduction and the valence bands near the bandgap but has more significant impact on the calculated bandgap values, which is due to the difference in the calculated lattice constants by different functionals.

Crystal structure relaxation and total energy calculations were performed using the Vienna ab initio simulation package (VASP).²⁸ Projector augmented wave (PAW)²⁹ pseudopotentials (treating 5s²5p⁶6s¹ of Cs, 4d¹⁰5s¹ of Ag, 5d¹⁰6s²6p³ of Bi, 3s²3p⁵ of Cl, and 4s²4p⁵ of Br as valence electrons) were used with the plane-wave basis cutoff energy set to 300 eV. The 6 × 6 × 6 Monkhorst-Pack k-point meshes were used for the Brillouin-zone integration of the primitive cell of Cs₂AgBiBr₆ and Cs₂AgBiCl₆. Adaptive k-points were used for secondary phases. We tested the convergence of the cutoff energy (up to 450 eV) and k-points (8 × 8 × 8 meshes) and found that the increased cutoff energy and large number of k-points produced negligible difference in the calculated results. We also performed test calculations with the 4p electrons of Ag as the valence electrons and found that the influence is negligible. Six different approximations to the exchange-correlation functionals, including LDA,¹⁸ two GGA functionals (PBE and PBE-sol),^19,24 and three van der Waals functionals (vdW-DF2, optB86b-vdW, and optB88-vdW),^25–27 were used as discussed as follows.

To quantify the error of the different exchange-correlation functionals, the statistical errors between theoretically calculated lattice constants and experimental data were computed, including the mean error (ME) and the mean absolute error (MAE), as well as the relative versions of these quantities, i.e., the mean relative error (MRE) and mean absolute relative error (MARE).

To evaluate the stability of Cs₂AgBiCl₆ and Cs₂AgBiBr₆, i.e., whether they can be synthesized in a certain chemical environment, we introduce chemical potentials of the component elements, which correspond to the chemical environment during synthesis. Under thermodynamic equilibrium conditions for crystal growth, the chemical potentials must meet several conditions in order to maintain the stable growth of the crystal and avoid competing phases. These conditions are shown below for Cs₂AgBiBr₆ as an example.

If the quaternary compound is stable during synthesis, the following thermodynamic equilibrium condition should be satisfied:

where μ_Cs, μ_Ag, μ_Bi, and μ_Br are the chemical potentials of Cs, Ag, Bi, and Br, respectively, relative to their elemental phases. ΔH_f(Cs₂AgBiBr₆) is the calculated formation energy of Cs₂AgBiBr₆. Note that all calculated formation energies in this work are given for per formula unit.

To exclude the coexistence of elemental phases, the following conditions should be satisfied:

To avoid the formation of the competing phases including CsBr₃, CsBr, AgBr, BiBr₃, BiBr, CsAgBr₂, CsAgBr₃, Cs₂AgBr₃, Cs₃Bi₂Br₉, and Ag₃(Bi₂Br₉)₇, the following constraints should be satisfied as well:

By solving the mathematical equation and inequation arrays, we can determine if there exists a certain chemical potential (μ_Cs, μ_Ag, μ_Bi, μ_Br) that can satisfy all these thermodynamic constraints [Eqs. (1)–(12)]. If there exists a range of μ_Cs, μ_Ag, μ_Bi, μ_Br satisfying the above constraints, Cs₂AgBiBr₆ is thermodynamically stable and single-phase crystals can be synthesized without the coexistence of any secondary phases in the environment corresponding to the chemical potentials within the range. If there is no stable region in the chemical potential space, Cs₂AgBiBr₆ is thermodynamically unstable and it will tend to phase-separate even if it is synthesized.

When we use the PBE functional to calculate the formation energies for all these compounds (for secondary phase, the lowest-energy ground-state structure according to the PBE calculations are used), we found that there does not exist a stable region in the chemical potential space that satisfies all the thermodynamic conditions, i.e., no chemical potentials can satisfy all the conditions in Eqs. (2)–(12) for Cs₂AgBiBr₆. The two-dimensional phase diagram with $μCs=−3.0eV$ is shown in Fig. 1(a) as an example. (Note that there are three independent chemical potentials for a quaternary compound. Here we choose μ_Cs, μ_Bi, and μ_Br as independent chemical potentials while μ_Ag can be related to the other three using Eq. (1). Equations (2)–(12) are evaluated in the three-dimensional chemical potential space to construct a phase diagram.) No other values of $μCs$ is found to lead to a stable region for Cs₂AgBiBr₆ in the phase diagram. This result indicates that Cs₂AgBiBr₆ is thermodynamically unstable with respect to other secondary phases.

In contrast to the predicted instability of Cs₂AgBiBr₆ by PBE calculations, large single crystals of Cs₂AgBiBr₆ have been synthesized experimentally.¹³ To understand the causes of the apparent contradiction, we performed additional calculations using LDA and the optB88-vdW functional that incorporates the van de Waals interaction. Interestingly, the use of both functionals leads to the prediction of a stable region for Cs₂AgBiBr₆ in the chemical potential space, as depicted in Figs. 1(b) and 1(c), indicating that the Cs₂AgBiBr₆ is thermodynamically stable. These results show that the prediction of thermodynamic stability by DFT calculations is sensitive to the exchange-correlation functionals.

PBE-sol, vdW-DF2, and optB86b-vdW functionals are further tested and the results are compared with those of LDA, PBE, and optB88-vdW calculations. To simplify the problem, we did not construct the three-dimensional phase diagram using each of these functionals. Instead, we calculated the decomposition energies of three possible phase-separation pathways of double-perovskite halides as reported in Ref. 18. These phase-separation pathways for Cs₂AgBiBr₆ are shown below. If the decomposition energy for any pathway is negative, no stable region for the double-perovskite halides in the phase diagram is possible,

The decomposition energies for the three pathways for both Cs₂AgBiCl₆ and Cs₂AgBiBr₆ are summarized in Table I. It can be seen that the decomposition energies of Cs₂AgBiBr₆ through the pathways 2 and 3 are negative according to the PBE calculations, indicating that Cs₂AgBiBr₆ is unstable and should tend to phase-separate through these two decomposition pathways spontaneously even if it is synthesized. This explains why there is no stable region for Cs₂AgBiBr₆ in the phase diagram in Fig. 1(a). By contrast, the calculated decomposed energies using the LDA and optB88-vdW functionals are positive, consistent with the calculated phase diagrams in Figs. 1(b) and 1(c).

The PBE calculations arrived at the wrong conclusion that Cs₂AgBiBr₆ is thermodynamically unstable, while the results obtained using all other functionals are consistent with the experiments. The failure of the PBE calculation can be traced back to the incorrect ground-state structure of the competing phase AgBr. When we calculated the formation energies for the competing phases, we considered all the possible crystal structures that were collected in the Materials Project database³⁰ and used the formation energy of the low-energy structure. According to our PBE calculations and also the PBE results collected in the Materials Project database, the lowest-energy (ground-state) structure of AgBr is the zinc blende structure (space group F-43m). In Fig. 2(a), we plotted the PBE-calculated relative total energies of the zinc blende and rock-salt structured AgBr as a function of the volume. The predicted lowest-energy structure (calculated ground-state) is the zinc blende structure, while the rock-salt structure (real ground-state) has a higher energy by 80 meV/f.u. at the equilibrium state. However, it is well known experimentally that the rock-salt structure (space group Fm-3m) is the ground-state structure for AgBr at room temperature.³¹ 

If we use the formation energy of the rock-salt structured AgBr (higher than that of the zinc blende AgBr in PBE calculations), we can find a stable region for Cs₂AgBiBr₆ in the chemical potential space, so it becomes stable and its decomposition energies also become positive, as shown by the values in the parentheses of Table I. In a previous study by Xiao et al.¹⁶ using the PBE functional, they also used the rock-salt structure and showed that Cs₂AgBiBr₆ is stable, in agreement with our calculations here. These results show that the instability of Cs₂AgBiBr₆ predicted by the PBE calculations results from the incorrect ground-state structure (zinc blende) of AgBr, which has lower energy than the rock-salt structure according to the PBE calculations.

By contrast, the LDA functional can predict the ground-state of AgBr correctly, i.e., the rock-salt structure has lower energy at equilibrium than the zinc blende structure, as shown in Fig. 2(b). Therefore, the rock-salt structure is used for predicting the stability, arriving at the correct conclusion that Cs₂AgBiBr₆ is stable as shown in Fig. 1(b) and Table I.

To verify whether the failure of PBE functional in predicting the ground-state structure is general for other binary halides, we calculated the relative total energy as a function of the volume for AgCl and CsCl using both PBE and LDA functionals. As shown in Figs. 2(c)–2(f), the PBE results contradict the LDA results for all three binary halides. For AgCl, the experimental ground-state is the rock-salt structure which is correctly predicted by the LDA calculation, but the PBE calculation shows that the zinc blende structure has lower energy at equilibrium. For CsCl, the experimental ground-state is obviously the CsCl structure which is correctly predicted by the LDA calculation, but the PBE calculation shows that the rock-salt structure has even lower energy.

Both previous calculations (e.g., the Materials Project database and Ref. 32) and our present calculations showed that the silver halides become more stable in the zinc blende structure than the rock-salt structure when the PBE functional is used. One possible explanation to this had been given in Ref. 32: The PBE calculations always prefer the more covalent zinc blende configurations with respect to the rock-salt structure. The distribution of electron densities for the rock-salt structure is more homogeneous than that for the zinc blende structure. As a consequence, the energy difference between zinc blende and rock-salt phases is enhanced in the PBE calculations (which favors the phase with larger inhomogeneity of the electron density), predicting the wrong ground-state structure for AgCl and AgBr.

These results indicate that the PBE functional is not accurate for studying the energetics and structural properties of binary metal halides and, therefore, is not appropriate for predicting the thermodynamic stability of double-perovskite halides, to which special attention should be paid in the future high-throughput Materials Genome study of metal halides.

Since the LDA and PBE functionals are inconsistent with each other in predicting low-energy structures of metal halides, we performed more test calculations using a series of different functionals, including LDA,¹⁸ two GGA forms (PBE and PBE-sol),^19,24 and three van der Waals functional forms (vdW-DF2, optB86b-vdW, and optB88-vdW).^25–27 We calculated total energies of CsCl, AgBr, and AgCl in two different crystal structures (the ground-state structure and an additional low-energy competing phase) as shown in Fig. 3. The calculated total energies for the competing phases are referenced to those of the experimental ground-state structures (denoted as real ground-state in Fig. 3), which are set to zero. The PBE functional predicts the incorrect ground-state structures for all three compounds, i.e., the calculated energies of the competing phases are lower than those of the real ground-state structures (negative relative energies). The PBE-sol and vdW-DF2 functionals predict the correct ground-states structure for both AgBr and AgCl but incorrect ones for CsCl, while the LDA, optB86b-vdW, and optB88-vdW functionals predict the correct ground-state structures for all three compounds.

We also compared the decomposition energies of Cs₂AgBiBr₆ and Cs₂AgBiCl₆ calculated using different functionals, as shown in Table I. Only the PBE calculations arrive at the incorrect conclusion that Cs₂AgBiBr₆ is unstable, while all other functionals arrive at the correct conclusion. Even the PBE-sol and vdW-DF2 functionals give correct results although they predict the incorrect ground-state structure for CsCl. Among the three more reliable functionals (LDA, optB86b-vdW, and optB88-vdW), the LDA functional gives relatively larger decomposition energies (stabilizing the quaternary compounds) than the two vdW functionals. Although there are differences in the absolute values, the trends in the LDA results are consistent with those in the optB86b-vdW and optB88-vdW results.

Using different exchange-correlation functionals in DFT calculations not only leads to qualitative differences in predicted thermodynamic stability of halides but also has a significant impact on the calculated lattice constants, which in turn affect the calculated bandgaps. The calculated and the experimental lattice constants of Cs₂AgBiCl₆, Cs₂AgBiBr₆, and their competing phases are summarized in Tables II and III. The mean error (ME), mean absolute error (MAE), mean relative error (MRE), and mean absolute relative error (MARE) between the calculated and the experimental values are also listed. Among the five functionals, the PBE and vdW-DF2 functionals overestimate lattice constant with positive MRE; e.g., their MRE is about 1.9% and 2.7%, respectively, for Cs₂AgBiBr₆ related compounds. By contrast, the LDA, PBE-sol, optB86b-vdW, and optB88-vdW all underestimate the lattice constants with negative MRE.

The errors (MARE) in the calculated lattice constants based on LDA calculations are 3.5% and 4.1% for Cs₂AgBiBr₆ and Cs₂AgBiCl₆, respectively. These errors are significant but may still be acceptable for many applications. The MAREs of PBE and PBE-sol calculations are significantly smaller than that of LDA while the vdW-DF2 calculations give somewhat larger errors than PBE calculations.

Using optB86b-vdW and optB88-vdW functionals consistently gives much smaller MAREs compared to other functionals. In particular, the lattice constants obtained by optB88-vdW calculations are the most accurate. Since the optB88-vdW functional gives excellent results on the ground-state structures, the stability trend, and the lattice constants for a range of binary, ternary, and quaternary metal halides, we conclude that the optB88-vdW functional should be reliable for studying energetic and structural properties of double-perovskite halides.

After comparing the accuracy of different functionals in calculating the formation energies and lattice constants, now we will analyze how the functionals influence the calculated electronic band structure.

The calculated band structures of Cs₂AgBiCl₆ and Cs₂AgBiBr₆ using different functionals are plotted in Figs. 4 and 5, respectively. It can be seen that using different functionals has negligible effects on the dispersion of the highest valence band and the lowest conduction bands. However, the bandgaps calculated by different functionals are significantly different; they vary from 1.42 eV to 1.90 eV for Cs₂AgBiCl₆ and from 0.93 eV to 1.36 eV for Cs₂AgBiBr₆. The PBE and vdW-DF2 functionals give the largest bandgaps, while the LDA functional always gives the smallest bandgaps for both Cs₂AgBiBr₆ and Cs₂AgBiCl₆. This can be ascribed to the overestimated lattice constants by PBE and vdW-DF2 and underestimated lattice constants by LDA. Previous studies had shown that the bandgaps of the Pb halides increase with their lattice constants due to the anti-bonding character of the valence band edge states (the antibonding state of the Pb 5s and I 5p hybridization which is weakened if the lattice constants are increased^5,40). The LDA functional predicts smaller lattice constants, so the bandgaps are smaller. The PBE and vdW-DF2 functionals predict larger lattice constants, so the bandgaps are larger. To demonstrate the influence of the lattice constant difference, we also calculated the band structures using different functionals for Cs₂AgBiCl₆ and Cs₂AgBiBr₆ with the same lattice constants (optimized by optB88-vdW). The band structures are shown in the supplementary material. As shown in Figs. S1 and S2 of the supplementary material, different functionals still have negligible effects on the dispersion of the highest valence band and the lowest conduction band; however, the bandgap difference calculated by different functionals becomes obviously smaller, i.e., the values vary from 1.56 eV to 1.82 eV for Cs₂AgBiCl₆ (compared to 1.42-1.90 eV when the lattice constants are different) and from 1.07 eV to 1.28 eV for Cs₂AgBiBr₆ (compared to 0.93-1.36 eV). Therefore, we can conclude that the differences in the lattice constants relaxed by different functionals play an important role in the large bandgap difference.

All the functionals predict indirect bandgaps for double-perovskite Cs₂AgBiCl₆ and Cs₂AgBiBr₆, which is consistent with the experimental observation.¹³ The values of the calculated bandgaps are smaller than the experimental values (2.77 eV for Cs₂AgBiCl₆ and 2.17 eV for Cs₂AgBiBr₆) because these functionals usually underestimate the bandgaps of semiconductors. Using hybrid functionals, which include a fraction of Fock exchange in the exchange functional, may increase the bandgaps.

Since the different functionals do not change the nature of the bandgap but influence the bandgap sizes significantly due to the lattice constant difference, we conclude that (i) the most accurate band structures and bandgaps should be calculated using the lattice constants relaxed by the optB88-vdW functional (which has the smallest error in the calculated lattice constant) and (ii) all the functionals can be used if only the band structure shape or bandgap nature is needed.

We found that the first-principles prediction of thermodynamic stability of double-perovskite halides is extremely sensitive to the approximations to the exchange-correlation functionals. The extensively used PBE functional gives an incorrect prediction that Cs₂AgBiBr₆ is thermodynamically unstable because the PBE calculation fails to predict the correct ground-state structure of AgBr, which is a competing phase of Cs₂AgBiBr₆. Such a failure of the PBE functional is found to exist also for AgCl and CsCl and, thus, may be general for halides. We systematically compared a series of functionals in predicting the thermodynamic stability of two double-perovskite halides, Cs₂AgBiBr₆ and Cs₂AgBiCl₆, and found that LDA, optB86b-vdW, and optB88-vdW can predict the ground-state structures of AgBr, AgCl, and CsCl correctly and also give correct trends consistently in the calculated decomposition energies. Furthermore, optB88-vdW gives a more accurate description of lattice constants than other functionals for a number of binary, ternary, and quaternary halides. Considering these two aspects, we conclude that the optB88-vdW functional is the most accurate and reliable for studying energetic and structural properties of double-perovskite halides and should be used in the future high-throughput calculation studies. Different functionals have little influence on the dispersion of the conduction and valence bands near the bandgaps, but the bandgap sizes calculated by different functionals differ significantly because the calculated lattice constants differ. Therefore, the most accurate band structures and bandgaps should be calculated using the lattice constants relaxed by the optB88-vdW functional.

See supplementary material for the band structures using different functionals for Cs₂AgBiCl₆ and Cs₂AgBiBr₆ with the same lattice constants (optimized by optB88-vdW).

This work was supported by the National Natural Science Foundation of China (NSFC) under Grant Nos. 61574059 and 61722402, National Key Research and Development Program of China (No. 2016YFB0700700), Shu-Guang Program (No. 15SG20), and CC of ECNU. Mao-Hua Du was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division.