Exchange-correlation functional challenges in modeling quaternary chalcogenides

The development of next-generation quaternary chalcogenides, such as $\mathrm{C}{\mathrm{u}}_{2}\mathrm{ZnSn}{\mathrm{S}}_{4}$ (CZTS) and $\mathrm{C}{\mathrm{u}}_{2}\mathrm{ZnGe}{\mathrm{S}}_{4}$ (CZGS), for solar energy and thermoelectric applications hinges upon both careful experimentation and accurate quantum mechanical modeling. To address the latter, many have turned to density functional theory (DFT), which offers several choices for the approximate treatment of electron exchange and correlation (XC). Popular XC functionals include the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) and the recently developed strongly constrained and appropriately normed (SCAN) meta-GGA. Extensions of DFT functionals, such as adding a Hubbard $U$ correction and introducing a fraction of the Fock exchange (hybrid functionals), have been used widely to model systems containing $3d$ metal ions. However, no studies yet have compared comprehensively PBE($+U$) and SCAN($+U$) in the quality of their predictions of the bulk and defect thermodynamics of quaternary chalcogenides, which play a critical role in device fabrication and performance. Hence, here we calculate the (i) 0 K formation energies of bulk Ge compounds and (ii) neutral defect formation energies including charge-balanced (e.g., $\mathrm{C}{\mathrm{u}}_{\mathrm{Zn}}+\mathrm{Z}{\mathrm{n}}_{\mathrm{Cu}}$) and charge-imbalanced (e.g., $\mathrm{C}{\mathrm{u}}_{\mathrm{Sn}}$) combinations of antisites and vacancies in CZTS and CZGS using the PBE, $\mathrm{PBE}+U$, SCAN, $\mathrm{SCAN}+U$, and the hybrid Heyd-Scuseria-Ernzerhof XC frameworks. We find that the formation energies of charge-imbalanced defects are more sensitive to the choice of the XC functional than those of charge-balanced defects, which can be explained by the differences in the extent of penalization of defect-generated delocalized electrons/holes by PBE, $\mathrm{PBE}+U$, SCAN, and $\mathrm{SCAN}+U$. Additionally, our results show that SCAN systematically underbinds Ge-containing compounds, thus highlighting the need for even further improvement of XC functionals. Based on our findings, we recommend the use of SCAN for modeling quaternary chalcogenides because its errors are systematic, and it has the firmest theoretical underpinning. Our work provides guidance for future modeling of quaternary chalcogenides.