Climbing Jacob's ladder: A density functional theory case study for Ag₂ZnSnSe₄ and Cu₂ZnSnSe₄

At ambient condition AZTSe and CZTSe crystallise in the kesterite crystal structure (space group $I\bar{4}$, no. 82), as shown in figure 1(a). The kesterite crystal structure can be most easily understood as originating from a binary II–VI zinc-blende crystal structure via pairwise cation substitutions, with the restriction that the bonding to adjacent cations has to fulfil the so-called octet rule, i.e. the valence shell of each atom comprises eight electrons. Starting from a II–VI zinc-blende material, doubling the unit cell along the c axis, and additionally replacing the group-II element by a pair of group-I and group-III elements, can lead to ternary I–III–VI₂ compounds of the chalcopyrite and the CuAu-like structure, respectively. A further replacement of the group-III elements by pairs of group-II and group-IV elements in the chalcopyrite structure leads to the I₂–II–IV–VI₄ kesterite crystal structure. However, performing the same replacement starting from the CuAu-like structure can lead to the I₂–II–IV–VI₄ stannite crystal structure (space group $I\bar{4}2{m}$, no. 121), which is shown in figure 1(c).

Zoom In Zoom Out Reset image size

To make matters worse, the ionic radii of Cu and Zn are nearly identical, giving rise to a high probability of cation exchange in the Cu–Zn planes perpendicular to the crystallographic c axis, located at z = 0.25 and z = 0.75. While partial exchange of Cu and Zn cations leads to an increase in the Cu$_{\textrm{Zn}}$ and Zn$_{\textrm{Cu}}$ defect density, a totally random distribution of Cu and Zn cations within the Cu–Zn planes leads to the introduction of additional symmetry elements and has therefore been named disordered kesterite structure, shown in figure 1(b). In terms of structural investigations this disordered kesterite structure crystallises in the same space group as the stannite crystal structure, and has hindered the correct crystal structure determination of possible solar absorber materials for some time.

In recent years, first-principles calculations based on DFT have become a very powerful tool for materials science investigations, not only due to the development of ever more accurate exchange and correlation functionals, but also due to the wide-spread availability of local, regional, and national high-performance computing (HPC) facilities, providing the necessary computational resources for investigations of increasing complexities. With the constant development of improved exchange and correlation functionals to be employed in first-principles investigations, for somebody not that familiar with the topic it became more and more difficult to judge the accuracy and reliability of the reported results. In order to allow for some guidance, the so-called Jacob's ladder can provide a first hint [16], as adopted in figure 2.

Zoom In Zoom Out Reset image size

Here, exchange and correlation functionals of similar capabilities are grouped together at the same rung of this ladder, while improved and more accurate ones are placed at a higher rung. Functionals of the local density approximation (LDA) can be found at the lowest rung and are based solely on the electron density, the more common examples include the parametrisations after Perdew and Zunger [17] and Perdew and Wang [18]. The next rung depicts the generalised gradient approximation (GGA) functionals, which additionally take into account the gradient of the electron density. Examples include the parametrisations after Perdew, Burke, and Ernzerhof (PBE) [19], the PBE parametrisation revised for solids (PBEsol) [20], and the parametrisation after Armiento and Mattsson (AM05) [21], respectively. Additional inclusion of the second derivative of the electron density (the Laplacian) and the orbital kinetic energy density leads to the rung of the so-called meta-GGAs, with the SCAN functional [22] being a prominent example. This SCAN functional satisfies all known possible exact constraints for the exact density functional, and has been claimed to match or improve on the accuracy of computationally more demanding hybrid functionals [23].

The mentioned hybrid functionals replace a pre-defined fraction of the underlying exchange and correlation energy by a Hartree–Fock exact exchange term. Prominent examples include the parametrisations after Heyd, Scuseria, and Ernzerhof [24], or the PBE0 functional [25]. Recent years witnessed a surge in investigations trying to improve hybrid functional investigations by adjusting the pre-defined fraction of Hartree–Fock exact exchange, ultimately being determined in a fully self-consistent manner [26, 27].

On the one hand, by climbing Jacob's ladder and employing ever more sophisticated exchange and correlation functionals, we expect to obtain improved accuracy from our first-principles calculations, as indicated by the right arrow in figure 2. On the other hand, this is accompanied by an increase in computational resources, as depicted by the left arrow in figure 2. However, with the more wide-spread availability of HPC facilities, hybrid functional investigations became more and more common in recent years and led to a better description of structural, electronic, and optical properties of a range of materials [28–31].

Apart from the calculation of electronic and optical properties by means of exchange and correlation functionals of four distinct rungs of Jacob's ladder as described above, we also performed quasiparticle calculations of those material properties based on the GW approximation introduced by Hedin [32]. While the GW approximation remains the most accurate method for quasiparticle calculations available, its numerical requirements and perturbative nature make approximations unavoidable [33–37].

From a numerical requirements' point of view several flavours of practical GW calculations have been established over the years. The most simple one is the single-shot $G_{0}W_{0}$ method, with G₀ being the noninteracting Green's function of the system and W₀ its screened Coulomb interaction. Extensions iterate the eigenvalues to self-consistency in the Green's function alone (GW₀), or also in the screened interaction W( GW). Methodical extensions see the inclusion of vertex corrections as well [38], however, this is beyond the scope of the present work.

From a perturbative nature's point of view, all mentioned flavours of GW calculations depend on the starting one-electron energies and orbitals, usually taken as Kohn-Sham eigenenergies and orbitals from preceding DFT calculations. This immediately implies that the accuracy of GW calculations will be influenced by the accuracy of the starting one-electron energies and orbitals, as will be discussed later in the results section.

All the calculations of the present work have been performed using the Vienna ab initio Simulations package (VASP, 5.4.4) [39–41] together with the projector-augmented wave (PAW) method [42, 43]. Structural relaxations employed the recommended PAW potentials supplied by VASP that provided 17, 17, 12, 14, and 6 valence electrons for the Cu, Ag, Zn, Sn, and Se atoms, respectively. The electronic band structures, the optical properties, and the many-body perturbation theory calculations based on the GW approximation made use of the PAW potentials recommended for GW calculations, and provided 19, 19, 20, 14, and 6 valence electrons for the Cu, Ag, Zn, Sn, and Se atoms instead.

In order to judge the accuracy and computational efficiency of several rungs of Jacob's ladder, we employed the following exchange-correlation potentials: for the LDA functional we chose the Perdew and Zunger parametrisation [17, 44], for the GGA functional we chose the Perdew, Burke, and Ernzerhof implementation revised for solids (PBEsol) [20], as meta-GGA functional we chose the relatively new SCAN functional [22], and finally as hybrid functional we chose the range-separated HSE06 functional [24], respectively.

Structural relaxations of the kesterite and stannite crystal structures have been performed for the standard primitive unit cells, containing one formula unit (f.u.) for both structural polymorphs. The ground state structures have been optimised by analysing the total energy curves, which have been obtained for several volumes around the experimentally known ground state volume. Keeping these volumes fixed, the respective lattice constants and all internal coordinates have been allowed to relax until the forces on all atoms were below 0.001 eV Å^–1. While the ground state volume V₀ has been obtained by a spline fit to the total energy curves, a subsequent analysis via Murnaghan's equation of state [45, 46]

$$\begin{equation} E(V) = E(V_{0}) + \frac{B_{0}V}{B_{0}^{\prime}} \left[ \frac{(V_{0}/V)^{B_{0}^{\prime}}}{B_{0}^{\prime}-1} +1 \right] - \frac{V_{0}B_{0}}{B_{0}^{\prime}-1}\end{equation} \tag{ 1 }$$

additionally yielded the bulk modulus B₀ and its pressure derivative $B_{0}^{\prime}$, respectively. A final relaxation of internal coordinates has been performed for the ground state volumes, yielding the overall ground states that have been further analysed using the AFLOW [47] and FINDSYM [48, 49] packages.

The obtained relaxed ground state structures served as a starting point for subsequent calculations of the electronic band structures and the real and imaginary parts of the dielectric functions, which have been obtained by summing over empty states using Fermi's Golden Rule, transition matrix elements, and applying a Kramers-Kronig transformation [50]. In order to ensure converged results the number of empty bands in the calculations of the optical properties have been increased by a factor of four. The final real and imaginary parts of the dielectric functions have been obtained by diagonalising the dielectric tensors for every energy point and averaging over the resulting main diagonal elements, as applied before to non-cubic oxide [51] and amorphous materials [52, 53].

Together with the other technical parameters, k-point grid of 6 × 6 × 6, cut-off energy of the plane-wave expansion of 500 eV, and a convergence criteria for the total energy of 10⁻⁶ eV, this ensured well-converged results. Due to the increased numerical demand the k-point grid has been reduced to 4 × 4 × 4 for the GW calculations.

The total energy curves for various employed exchange and correlation functionals, normalised to one functional unit (f.u.) and rescaled to zero energy, are shown in figure 3 for the kesterite crystal structures of AZTSe (a) and CZTSe (b), respectively. The dashed vertical lines indicate the experimental ground state volumes taken from a combined x-ray and neutron powder diffraction study of Gurieva et al[ 2]. The respective total energy curves for the stannite crystal structures of AZTSe and CZTSe can be found in the appendix (figure A1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The obtained ground state structural properties of the kesterite and stannite crystal structures of AZTSe and CZTSe are given in tables 1 and 2 for various exchange and correlation functionals in comparison to experimental data, respectively. For both, kesterite-type AZTSe and CZTSe, the LDA calculated ground state volume is severely underestimated by approximately 4%, as shown in figures 3 and 4. Deviations from the experimental ground state volumes of about 1% are obtained for the PBEsol and the SCAN functionals, with a much better agreement of the SCAN functional in case of kesterite-type CZTSe. The hybrid functional HSE06 overestimates the experimental ground state volumes by 2%–3%, respectively.

Zoom In Zoom Out Reset image size

Table 1. Ground state structural properties (ground state volume V₀, lattice constants a₀ and c₀, fractional coordinates x, y, and z of the anion Wyckoff position 8g, and bulk modulus B₀ and its pressure derivative $B_{0}^{\prime}$) obtained by various exchange and correlation functionals in comparison to experimental data for the kesterite crystal structures of Ag₂ZnSnSe₄ and Cu₂ZnSnSe₄.

Ag2ZnSnSe4	V0( Å3)	a0( Å)	c0( Å)	x	y	z	B0( GPa)	$B_{0}^{\prime}$
LDA	198.912	5.916	11.365	0.262	0.241	0.131	62.8	5.0
PBEsol	205.756	5.987	11.479	0.264	0.242	0.130	56.1	4.9
SCAN	209.448	6.072	11.362	0.269	0.243	0.130	53.5	5.0
HSE06	213.628	6.096	11.496	0.272	0.246	0.128	49.3	5.0
exp. [2]	207.262	6.046	11.340
exp. [11]	205.866	6.036	11.301	0.251	0.231	0.128
Cu2ZnSnSe4
LDA	176.011	5.605	11.207	0.228	0.239	0.131	72.9	4.8
PBEsol	181.736	5.665	11.327	0.229	0.239	0.131	66.1	5.0
SCAN	183.887	5.695	11.340	0.232	0.240	0.130	64.3	5.2
HSE06	187.615	5.732	11.422	0.239	0.244	0.129	59.6	5.0
exp. [1]	183.976	5.695	11.345
exp. [11]	183.702	5.692	11.340	0.251	0.240	0.128
exp. [2]	183.880	5.693	11.347

Table 2. Ground state structural properties (ground state volume V₀, lattice constants a₀ and c₀, fractional coordinates x and z of the anion Wyckoff position 8i, and bulk modulus B₀ and its pressure derivative $B_{0}^{\prime}$) obtained by various exchange and correlation functionals for the stannite crystal structures of Ag₂ZnSnSe₄ and Cu₂ZnSnSe₄.

Ag2ZnSnSe4	V0( Å3)	a0( Å)	c0( Å)	x	z	B0( GPa)	$B_{0}^{\prime}$
LDA	199.001	5.724	12.149	0.759	0.879	62.4	4.6
PBEsol	205.720	5.776	12.334	0.758	0.880	55.7	5.0
SCAN	209.603	5.768	12.599	0.758	0.882	53.6	3.3
HSE06	213.282	5.800	12.681	0.754	0.884	49.6	4.9
Cu2ZnSnSe4
LDA	176.051	5.607	11.198	0.761	0.863	72.8	4.7
PBEsol	181.793	5.667	11.320	0.760	0.864	66.0	4.9
SCAN	183.838	5.677	11.407	0.760	0.865	64.7	5.0
HSE06	187.463	5.706	11.514	0.756	0.868	60.1	5.0

The obtained ground state volumes increase with the increasing sophistication of the employed exchange and correlation functionals (LDA, PBEsol, SCAN, and HSE06). This trend is equally observed for the lattice parameters a and c for both materials, AZTSe and CZTSe, and for both crystal structures, given for kesterite in table 1 and for stannite in table 2, respectively. The bulk modulus B₀ shows the same decreasing trend for both materials and crystal structures, and its pressure derivative is mostly around 5.

The fractional coordinates of the anion Wyckoff positions show a complementary trend, i.e. increasing (decreasing) for the kesterite-type 8g fractional coordinates x and y( z), and decreasing (increasing) for the stannite-type 8i fractional coordinates x( z), respectively. The performance of various exchange and correlation functionals in describing fractional coordinates will become important in comparison to precise experimental local structure investigations, and has been shown to substantially influence the electronic band gaps in CZTS and CZTSe [54].

Based on the obtained ground state structures we calculated the electronic band structures as outlined in section 2.3. The electronic band structures are shown in figure 5 for the kesterite crystal structures of AZTSe (a) and CZTSe (b), respectively. The respective electronic band structures for the stannite crystal structures of AZTSe and CZTSe can be found in the appendix (figure A2). In all electronic band structure figures, green and red lines depict the valence and conduction bands calculated with the hybrid HSE06 functional, whereas the shaded grey background indicates the SCAN calculated valence and conduction bands. The dotted black lines are the results of many-body perturbation theory calculations by means of the single-shot $G_{0}W_{0}$ approximation with the HSE06 eigenvalues and orbitals as starting points.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a first approximation to the electronic band gaps of the materials, we look at the Kohn-Sham eigenvalue differences between the valence and conduction bands for the LDA, PBEsol, SCAN, and HSE06 functionals, and the quasiparticle energy differences in case of $G_{0}W_{0}$ calculations, respectively. For kesterite-type AZTSe the electronic band gaps amount to 0.16, 0.12, 0.30, and 1.04 eV for the LDA, PBEsol, SCAN, and hybrid HSE06 functionals. The closest agreement to the experimental band gap of 1.31 eV [2] is obtained for the hybrid HSE06 calculation. Similarly, for kesterite-type CZTSe the electronic band gaps amount to 0.01, 0.02, and 0.03 eV, i.e. close to a metallic solution, and 0.87 eV for the LDA, PBEsol, SCAN, and hybrid HSE06 functionals. Again, the hybrid HSE06 calculation yields closest agreement with the experimental band gap ranging between 0.94 eV [2] and 1.0 eV [1]. The closer agreement of the calculated band gap energies with the experimental values in case of CZTSe is possibly due to the better agreement of the obtained fractional coordinates of the anion Wyckoff position with respect to the experimental values, as given in table 1.

From the band gap values and the electronic band structures in Figure 5 it becomes apparent, that the simpler non-hybrid exchange and correlation functionals show the known band gap underestimation with respect to the experimental values. This severe underestimation of the electronic band gaps obtained by non-hybrid functionals is shown in figure 6.

Zoom In Zoom Out Reset image size

One important objective of the present work is to benchmark different computational approaches against the agreement with experimental structural and electronic properties. As for the structural properties, overall best agreement with available experimental data has been obtained by the SCAN functional, and the PBEsol functional being very close. The electronic properties, however, show best agreement for the hybrid HSE06 functional.

While the present investigation deals with relatively small unit cells, where even a full structural relaxation employing the hybrid HSE06 functional is feasible, for many other possible investigations structural relaxations based on hybrid functional become prohibitively expensive in terms of computational resources. This includes all types of investigations requiring larger supercells, e.g. formation energies of defects and defect complexes or the calculation of band offsets in modelled device heterostructures. For these types of investigations it is indispensable to estimate the deviations introduced by less accurate exchange and correlation functionals.

Over the last years it became more and more popular to perform structural relaxations of larger unit cells employing computationally cheaper functionals, and obtain the electronic properties by means of hybrid functional calculations without further relaxing the structure, so-called single-shot or one-shot hybrid calculations. The HSE06 calculated electronic band gaps based on the SCAN optimised structures amount to 0.97 and 0.74 eV for AZTSe and CZTSe, being somewhat lower that the HSE06 values obtained for the HSE06 optimised structures of 1.04 and 0.87 eV, respectively. Compared to the full hybrid HSE06 investigation, this approach yields the second best agreement with respect to the experimental band gaps, however, for larger unit cells it is presently the only approach available for restricted computational resources.

In addition, we performed quasiparticle calculations based on the single-shot $G_{0}W_{0}$ method. For kesterite-type AZTSe and depending on the starting eigenvalues and orbitals, the $G_{0}W_{0}$ calculations yield a quasiparticle gap of 1.48, 1.47, and 1.56 eV for the plain SCAN functional, the HSE06 calculations on top of the SCAN optimised structures, and the hybrid HSE06 functional, respectively. The respective $G_{0}W_{0}$ calculations for kesterite-type CZTSe yield quasiparticle gaps of 1.43, 0.97, and 1.12 eV, with starting eigenvalues and orbitals from the plain SCAN functional, the HSE06 calculations on top of the SCAN optimised structures, and the hybrid HSE06 functional, respectively. It can be seen, that in all cases, the $G_{0}W_{0}$ quasiparticle gaps are overestimated with respect to the experimental band gaps; however, the agreement is best for the starting eigenvalues and orbitals taken from a single-shot hybrid HSE06 calculation on top of the SCAN optimised structures (figure 6). Various test calculations to iterate the Green's functions to self-consistency, i.e. GW₀ calculations, showed only negligible influence on the quasiparticle gaps.

In terms of a direct comparison between different theoretical approaches, the $G_{0}W_{0}$ quasiparticle gaps are enlarged by approximately 0.5 eV (0.25 eV) for kesterite-type AZTSe (CZTSe) for starting eigenvalues and orbitals taken from both, the HSE06 calculations on top of the SCAN optimised structures, and the hybrid HSE06 functional, respectively. The same trend is observed for the stannite-type structures. The total overestimation of the quasiparticle gap of 1.43 eV with starting eigenvalues and orbitals from a SCAN calculation for kesterite-type CZTSe highlights the perturbative character of the GW method in general, and that the near metallic solution of the SCAN calculation is inadequate to provide starting conditions for subsequent quasiparticle calculations.

Based on the overall best agreement with the experimental band gaps for the hybrid HSE06 calculations, we also calculated the optical properties. The real (red) and imaginary (green) parts of the dielectric function are shown in figure 7 for the kesterite crystal structures of AZTSe (a) and CZTSe (b), respectively. In case of CZTSe, the dashed lines present experimental results of Leon et al[ 55], obtained via spectroscopic ellipsometry on bulk crystals. Although the overall structure of the theoretical results is much more pronounced, the important peak positions agree well with each other. $G_{0}W_{0}$ calculations of the dielectric functions basically show the same peak structure, and are only shifted to reflect the change in the band gap energies. The respective dielectric functions for the stannite crystal structures of AZTSe and CZTSe can be found in the appendix (figure A3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size